# Should I use premium Diesel? Result: No

[This article was first published on

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

A while ago I had a post: ‘Should I use premium Diesel? Setup. Since that time the data has been acquired. This post describes the results.**Wiekvoet**, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

### Data

Data is registered by me in 2014 and 2015. 2014 has standard Diesel, while 2015 has premium. Both are from the same months of the year, so should have approximately the same weather and traffic.### Results

The figure shows the main result: While the 2015 has some better results, it is not by much. It is not in the ballpark of 5% which was originally defined.Results from distribution fitting and t.test confirm the visual observation.

#### fitdistr

Standardmean sd

3.59517971 0.19314598

(0.04828649) (0.03414371)

Premium

mean sd

3.54015573 0.23905047

(0.05976262) (0.04225855)

#### t.test

Welch Two Sample t-testdata: usage by kind

t = -0.6934, df = 28.732, p-value = 0.4936

alternative hypothesis: true difference in means is not equal to 0

95 percent confidence interval:

-0.2173822 0.1073343

sample estimates:

mean in group Premium mean in group Standard

3.540156 3.595180

#### Posterior density for premium Diesel

Previously I defined a prior for the analysis of the data. This was a mixed distribution from three parts. Like Standard fuel, 5% better and something else in the neighborhood of 3.6. Given the simplicity of the model and the non-standard prior I used MCMCpack’s MCMCmetrop1R function. This allowed me to write it all as R code.The result is a density for mean usage quite close to 3.6.

Iterations = 501:20500

Thinning interval = 1

Number of chains = 1

Sample size per chain = 20000

1. Empirical mean and standard deviation for each variable,

plus standard error of the mean:

Mean SD Naive SE Time-series SE

[1,] 3.5409 0.06413 0.0004535 0.001722

[2,] 0.2699 0.05353 0.0003785 0.001448

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%

var1 3.4147 3.494 3.5465 3.5884 3.6523

var2 0.1873 0.232 0.2629 0.2998 0.3953

### Conclusion

While premium Diesel had a few occasions with remarkable fuel usage, it did not result in a relevant decrease in fuel usage overall. In addition, while the advised prices for standard and premium are 5% apart, the achievable price difference is 10%. The real price fighters do not do premium, hence the pressure to discount is less. This 10% price increase is more than I initially expected. Hence the data do not show benefit for using premium Diesel.

### Data and code

library(MASS)

library(ggplot2)

library(MCMCpack)

r1 <- read.table(text='l km

28.25 710.6

22.93 690.4

28.51 760.5

23.22 697.9

31.52 871.2

24.68 689.6

30.85 826.9

23.04 699

29.96 845.3

30.16 894.7

25.71 696

23.6 669.8

28.57 739

27.23 727.4

18.31 499.9

24.28 689.5′,header=TRUE)

r1$usage=100*r1$l/r1$km

r1$Fuel=’Standard’

r2 <- read.table(text='l km

23.97 696

25.33 693.5

24.19 698.5

23.27 689.6

25.78 707.8

15.84 510.1

32.23 860.2

23.29 674.6

25.33 678.5

26.44 686

36.43 913.2

26.46 850.1

24.38 678.8

34.85 1001.3

28.34 858.5

25.85 698.8′

,header=TRUE)

r2$usage=100*r2$l/r2$km

r2$Fuel=’Premium’

rr <- rbind(r1,r2)

rr$Fuel <- factor(rr$Fuel)

ggplot(rr,aes(Fuel,usage))+geom_violin()+scale_y_continuous(‘usage (l/100 km)’)

fitdistr(r1$usage,’normal’)

fitdistr(r2$usage,’normal’)

t.test(usage ~ kind, data=rr)

priormeandens <- function(usage) {

dens <- (dnorm(usage,3.6,.05)+

dnorm(usage,3.6/1.05,.05)+

dnorm(usage,(3.6+3.6/1.05)/2,.15))/3

log(dens)

}

mydens <- function(theta,x) {

musage=theta[1]

sdusage=theta[2]

datdens <- sum(dnorm(x,musage,sdusage,log=TRUE))

priormean <- priormeandens(musage)

priorsd <- dunif(sdusage,0,1,log=TRUE)

datdens+priormean+priorsd

}

mcmc1 <- MCMCmetrop1R(mydens,

theta.init=c(3.8,.2),

seed=as.integer((Sys.time())),

x=r2$usage)

plot(mcmc1)

summary(mcmc1)

To

**leave a comment**for the author, please follow the link and comment on their blog:**Wiekvoet**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.

Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.