Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
0 parents
commit 2609718
Showing
390 changed files
with
567,792 additions
and
0 deletions.
There are no files selected for viewing
Binary file not shown.
Empty file.
Large diffs are not rendered by default.
Oops, something went wrong.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,37 @@ | ||
| ########## | ||
| # Compare Normal approximation of Beta posterior to the true posterior (MathExercise 3) | ||
| # | ||
| # Author: Per Siden | ||
| # Date: 2018-04-11 | ||
| ########## | ||
|
|
||
| # 3c | ||
| alpha = 1 | ||
| beta = 1 | ||
| n = 6 | ||
| s = 1 | ||
| f = n-s | ||
| normalMean = (alpha + s - 1)/(alpha+beta+n-2) | ||
| normalVar = (alpha + s - 1)*(beta + f - 1)/(alpha + beta + n - 2)^3 | ||
|
|
||
| xgrid = seq(0,1,.001) | ||
| plot(xgrid,dbeta(xgrid,alpha+s,beta+f),type="l") | ||
| lines(xgrid,dnorm(xgrid,normalMean,sqrt(normalVar)),col=2) | ||
| legend("topright", box.lty = 1, legend = c("Beta posterior","Normal approx."), | ||
| col = c("black",'red'), lwd = 2) | ||
|
|
||
| # 3d | ||
| alpha = 1 | ||
| beta = 1 | ||
| n = 12 | ||
| s = 2 | ||
| f = n-s | ||
| normalMean = (alpha + s - 1)/(alpha+beta+n-2) | ||
| normalVar = (alpha + s - 1)*(beta + f - 1)/(alpha + beta + n - 2)^3 | ||
|
|
||
| xgrid = seq(0,1,.001) | ||
| plot(xgrid,dbeta(xgrid,alpha+s,beta+f),type="l") | ||
| lines(xgrid,dnorm(xgrid,normalMean,sqrt(normalVar)),col=2) | ||
| legend("topright", box.lty = 1, legend = c("Beta posterior","Normal approx."), | ||
| col = c("black",'red'), lwd = 2) | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,206 @@ | ||
| logWage age | ||
| 11.1563 21.0000 | ||
| 12.8131 22.0000 | ||
| 13.0960 22.0000 | ||
| 11.6952 22.0000 | ||
| 11.5327 22.0000 | ||
| 12.7657 22.0000 | ||
| 12.5879 22.0000 | ||
| 11.9829 22.0000 | ||
| 13.4588 22.0000 | ||
| 12.2061 23.0000 | ||
| 12.0436 23.0000 | ||
| 11.9250 23.0000 | ||
| 13.7001 23.0000 | ||
| 12.7939 23.0000 | ||
| 13.3471 23.0000 | ||
| 13.0562 23.0000 | ||
| 12.0668 23.0000 | ||
| 13.2879 23.0000 | ||
| 13.6808 24.0000 | ||
| 12.8992 24.0000 | ||
| 12.2061 24.0000 | ||
| 13.5008 24.0000 | ||
| 13.5924 24.0000 | ||
| 13.6774 24.0000 | ||
| 13.4298 24.0000 | ||
| 13.1993 24.0000 | ||
| 13.0368 24.0000 | ||
| 11.8277 24.0000 | ||
| 13.3535 24.0000 | ||
| 12.5981 24.0000 | ||
| 13.4786 25.0000 | ||
| 13.2534 25.0000 | ||
| 13.4284 25.0000 | ||
| 13.7102 25.0000 | ||
| 13.6635 25.0000 | ||
| 13.5924 25.0000 | ||
| 13.0815 25.0000 | ||
| 13.0836 25.0000 | ||
| 13.3519 25.0000 | ||
| 13.2249 25.0000 | ||
| 13.4313 26.0000 | ||
| 13.3342 26.0000 | ||
| 13.1993 26.0000 | ||
| 13.3847 26.0000 | ||
| 13.3047 26.0000 | ||
| 14.0233 26.0000 | ||
| 13.4800 26.0000 | ||
| 13.8353 27.0000 | ||
| 13.9457 27.0000 | ||
| 13.3661 27.0000 | ||
| 13.5278 27.0000 | ||
| 13.3047 27.0000 | ||
| 13.7102 27.0000 | ||
| 13.7102 27.0000 | ||
| 13.5411 28.0000 | ||
| 13.7212 28.0000 | ||
| 13.5658 28.0000 | ||
| 13.3359 28.0000 | ||
| 13.8869 29.0000 | ||
| 13.9810 29.0000 | ||
| 13.5158 30.0000 | ||
| 13.6292 30.0000 | ||
| 13.3047 30.0000 | ||
| 13.4150 30.0000 | ||
| 13.6797 30.0000 | ||
| 13.8643 31.0000 | ||
| 13.4588 31.0000 | ||
| 13.5899 31.0000 | ||
| 13.7113 31.0000 | ||
| 13.4588 31.0000 | ||
| 13.5411 32.0000 | ||
| 13.3047 32.0000 | ||
| 13.9978 32.0000 | ||
| 13.4603 32.0000 | ||
| 13.8925 32.0000 | ||
| 13.8165 32.0000 | ||
| 13.8738 32.0000 | ||
| 13.8451 33.0000 | ||
| 13.5924 33.0000 | ||
| 13.4488 33.0000 | ||
| 13.9622 33.0000 | ||
| 13.6956 33.0000 | ||
| 13.0170 34.0000 | ||
| 13.9377 34.0000 | ||
| 14.0103 34.0000 | ||
| 13.3047 34.0000 | ||
| 14.1527 35.0000 | ||
| 13.8155 35.0000 | ||
| 13.9553 35.0000 | ||
| 13.9978 35.0000 | ||
| 13.7788 35.0000 | ||
| 13.8304 35.0000 | ||
| 13.7642 36.0000 | ||
| 13.4150 36.0000 | ||
| 13.9810 36.0000 | ||
| 13.9987 36.0000 | ||
| 13.6831 37.0000 | ||
| 13.0836 37.0000 | ||
| 14.1871 37.0000 | ||
| 14.1520 38.0000 | ||
| 13.5949 38.0000 | ||
| 13.4588 38.0000 | ||
| 13.8633 38.0000 | ||
| 14.1193 38.0000 | ||
| 14.0330 38.0000 | ||
| 13.2067 38.0000 | ||
| 13.6624 38.0000 | ||
| 13.8509 39.0000 | ||
| 13.4617 39.0000 | ||
| 14.0012 39.0000 | ||
| 13.5411 39.0000 | ||
| 14.1520 40.0000 | ||
| 13.6292 40.0000 | ||
| 13.7102 40.0000 | ||
| 14.2210 40.0000 | ||
| 13.3407 40.0000 | ||
| 13.8155 40.0000 | ||
| 13.1863 40.0000 | ||
| 13.8822 41.0000 | ||
| 13.8155 41.0000 | ||
| 14.2698 41.0000 | ||
| 12.3194 41.0000 | ||
| 13.9987 41.0000 | ||
| 12.3884 41.0000 | ||
| 13.7451 42.0000 | ||
| 14.1520 42.0000 | ||
| 13.9044 42.0000 | ||
| 13.7589 42.0000 | ||
| 13.7504 42.0000 | ||
| 13.5076 43.0000 | ||
| 13.7201 43.0000 | ||
| 13.4588 44.0000 | ||
| 14.6622 44.0000 | ||
| 12.5062 44.0000 | ||
| 12.9408 44.0000 | ||
| 13.9518 45.0000 | ||
| 13.6268 45.0000 | ||
| 13.8075 45.0000 | ||
| 13.3816 45.0000 | ||
| 13.8165 46.0000 | ||
| 13.9108 46.0000 | ||
| 13.8155 47.0000 | ||
| 13.7621 47.0000 | ||
| 13.7321 47.0000 | ||
| 13.8155 47.0000 | ||
| 13.4660 47.0000 | ||
| 11.4076 47.0000 | ||
| 13.9978 47.0000 | ||
| 14.0585 48.0000 | ||
| 12.8992 48.0000 | ||
| 14.0306 48.0000 | ||
| 13.3064 48.0000 | ||
| 13.2534 48.0000 | ||
| 13.9987 49.0000 | ||
| 13.9404 49.0000 | ||
| 13.9421 50.0000 | ||
| 13.4588 50.0000 | ||
| 13.7012 50.0000 | ||
| 14.1052 50.0000 | ||
| 14.6340 50.0000 | ||
| 13.4800 50.0000 | ||
| 14.4991 50.0000 | ||
| 13.6060 51.0000 | ||
| 13.9492 52.0000 | ||
| 13.7736 52.0000 | ||
| 13.8314 52.0000 | ||
| 14.0779 52.0000 | ||
| 14.2150 52.0000 | ||
| 13.6412 52.0000 | ||
| 11.8056 52.0000 | ||
| 14.2544 53.0000 | ||
| 15.0582 53.0000 | ||
| 13.2067 53.0000 | ||
| 13.6530 53.0000 | ||
| 13.0562 54.0000 | ||
| 13.8614 54.0000 | ||
| 14.2615 54.0000 | ||
| 14.7318 54.0000 | ||
| 13.4603 55.0000 | ||
| 13.6219 55.0000 | ||
| 13.2879 55.0000 | ||
| 13.5144 55.0000 | ||
| 14.4226 55.0000 | ||
| 13.6542 56.0000 | ||
| 13.7963 56.0000 | ||
| 12.2737 57.0000 | ||
| 14.3538 58.0000 | ||
| 13.3692 58.0000 | ||
| 13.3359 58.0000 | ||
| 13.3229 58.0000 | ||
| 12.8992 59.0000 | ||
| 13.9631 59.0000 | ||
| 13.8095 59.0000 | ||
| 14.5087 60.0000 | ||
| 11.5991 60.0000 | ||
| 13.4716 60.0000 | ||
| 13.5278 60.0000 | ||
| 12.5879 61.0000 | ||
| 13.3391 62.0000 | ||
| 12.6947 62.0000 | ||
| 13.9682 63.0000 | ||
| 12.0782 63.0000 | ||
| 14.0916 63.0000 | ||
| 13.7102 64.0000 | ||
| 12.2549 65.0000 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,43 @@ | ||
| ################################################################################################# | ||
| ################# Example of conjugate prior inference of the multinomial model ################ | ||
| ################################################################################################# | ||
|
|
||
| ####### Defining a function that simulates from a Dirichlet distribution | ||
| SimDirichlet <- function(nIter,param){ | ||
| nCat <- length(param) | ||
| thetaDraws <- matrix(0,nIter,nCat) # Matrix where the posterior draws are stored | ||
| for (j in 1:nCat){ | ||
| thetaDraws[,j] <- rgamma(nIter,param[j],1) | ||
| } | ||
| for (i in 1:nIter){ | ||
| thetaDraws[i,] = thetaDraws[i,]/sum(thetaDraws[i,]) # Diving every column of ThetaDraws by the sum of the elements in that column. | ||
| } | ||
| return(thetaDraws) | ||
| } | ||
|
|
||
| ########### Setting up data and prior ################# | ||
| y <- c(36,87,77) # Data | ||
| alpha <- c(1,1,1) # Dirichlet prior hyperparameters | ||
| nIter <- 10000 # Number of posterior draws | ||
|
|
||
| ########### Posterior sampling from Dirichlet ################# | ||
| thetaDraws <- SimDirichlet(nIter,y + alpha) | ||
|
|
||
| ################ Computing Summary statistics from the posterior sample ################### | ||
| mean(thetaDraws[,1]) | ||
| mean(thetaDraws[,2]) | ||
| mean(thetaDraws[,3]) | ||
|
|
||
| sqrt(var(thetaDraws[,1])) | ||
| sqrt(var(thetaDraws[,2])) | ||
| sqrt(var(thetaDraws[,3])) | ||
|
|
||
| sum(thetaDraws[,2]>thetaDraws[,3])/nIter # p(theta2>theta3|Data) | ||
|
|
||
| # Plots histograms of the posterior draws | ||
| plot.new() # Opens a new graphical window | ||
| par(mfrow = c(2,2)) # Splits the graphical window in four parts (2-by-2 structure) | ||
| hist(thetaDraws[,1],25) # Plots the histogram of theta[,1] in the upper left subgraph | ||
| hist(thetaDraws[,2],25) | ||
| hist(thetaDraws[,3],25) | ||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,26 @@ | ||
| y <- c(180,230,62,41) # The cell phone survey data (K=4) | ||
| alpha <- c(15,15,10,10) # Dirichlet prior hyperparameters | ||
| nIter <- 1000 # Number of posterior draws | ||
|
|
||
| # Defining a function that simulates from a Dirichlet distribution | ||
| SimDirichlet <- function(nIter, param){ | ||
| nCat <- length(param) | ||
| thetaDraws <- as.data.frame(matrix(NA, nIter, nCat)) # Storage. | ||
| for (j in 1:nCat){ | ||
| thetaDraws[,j] <- rgamma(nIter,param[j],1) | ||
| } | ||
| for (i in 1:nIter){ | ||
| thetaDraws[i,] = thetaDraws[i,]/sum(thetaDraws[i,]) | ||
| } | ||
| return(thetaDraws) | ||
| } | ||
| # Posterior sampling from Dirichlet posterior | ||
| thetaDraws <- SimDirichlet(nIter,y + alpha) | ||
|
|
||
| # Posterior mean and standard deviation of Androids share (in %) | ||
| message(mean(100*thetaDraws[,2])) | ||
| message(sd(100*thetaDraws[,2])) | ||
|
|
||
| # Computing the posterior probability that Android is the largest | ||
| PrAndroidLargest <- sum(thetaDraws[,2]>apply(thetaDraws[,c(1,3,4)],1,max))/nIter | ||
| message(paste('Pr(Android has the largest market share) = ', PrAndroidLargest)) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,16 @@ | ||
| # Bivariate density estimator | ||
| #install.packages("ks") # bivariate kernel density estimates | ||
| par(mfrow=c(1,2)) | ||
| library(ks) | ||
| Hpi1 <- Hscv(x = directDraws) | ||
| fEst <- kde(x = directDraws, H=Hpi1) | ||
| plot(fEst) | ||
| evalPoints <- t(rbind(fEst$eval.points[[1]],fEst$eval.points[[1]])) # Finding the points where the density estimate was evaluated in the plot | ||
|
|
||
| # Trick to compute the true density | ||
| bivariate <- function(x,y){ | ||
| return (dmvnorm(c(x,y), mean = mu, sigma = Sigma)) | ||
| } | ||
| trueDensity <- outer(evalPoints[,1],evalPoints[,2],bivariate) | ||
| contour(evalPoints[,1], evalPoints[,2], trueDensity, theta = 55, phi = 30, r = 40, d = 0.1, expand = 0.5, | ||
| ltheta = 90, lphi = 180, shade = 0.4, ticktype = "detailed", nticks=5) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,20 @@ | ||
| # Eight schools example in Rstan Vignette (https://cran.r-project.org/web/packages/rstan/vignettes/rstan.html) | ||
|
|
||
| data { | ||
| int<lower=0> J; // number of schools | ||
| real y[J]; // estimated treatment effects | ||
| real<lower=0> sigma[J]; // s.e. of effect estimates | ||
| } | ||
| parameters { | ||
| real mu; | ||
| real<lower=0> tau; | ||
| vector[J] eta; | ||
| } | ||
| transformed parameters { | ||
| vector[J] theta; | ||
| theta = mu + tau * eta; | ||
| } | ||
| model { | ||
| target += normal_lpdf(eta | 0, 1); | ||
| target += normal_lpdf(y | theta, sigma); | ||
| } |
Oops, something went wrong.