library('rstan')
library('fields')
library('parallel')

stan_model_file = 'monthly_precip_amount_gamma.stan'
iterations = 1500
warmup = 500
chains = 1
cores = 1#min(detectCores(), 4)

data_url = 'http://cires1.colorado.edu/~aslater/CVEN_6833/colo_monthly_precip_01.dat'
df = na.omit(read.table(data_url,na.strings='-999.999'))[1:100,]
names(df) <- c('lat','lon','elev','precip')


# first column is intercept
# scale the other covariates to reduce roundoff error
covars = cbind(1,scale(df[,c('lon','lat','elev')]))
y = as.vector(scale(df$precip,center=F))

S = length(y)
P = ncol(covars) 

# distance matrix in kilometers
dist = rdist.earth(df[,c('lon','lat')],miles=FALSE)

stan_data = list(S=S,P=P,y=y,covars=covars)

# initial guesses
initf = function()list(
    sill = 10,
    range = max(dist),
    nugget = 10,
    betas = coefficients(glm(y~covars+0, family=Gamma(link="log")))
    )

stanfit = stan(stan_model_file, data=stan_data, chains=chains, 
    iter=iterations, warmup=warmup, cores=cores, init=initf,
    refresh=1)
	
	
	### get the output stats - mean, median standard deviation from the posterior
### of the paramters
print(stanfit,digits=2)

# look at the results - in terms of convergence
traceplot(stanfit)

### Plot the posterior
dev.new()
par(mfcol=c(1,2))
hist(extract(stanfit)$betas[,1])
hist(extract(stanfit)$shape)
hist(extract(stanfit)$nugget)

### etc.



### use shinystan
launch_shinystan(stanfit)

#####  Posterior distribution