#!/bin/python# This code is based upon the exerice 4 from the IMM DTU course 02457# written by Lars Kai Hansen & Karam Sidarosfromnumpyimport*importnumpyasnpfromscipyimport*importpypr.clustering.gmmasgmmimportpypr.preprocessing.lag_matrixaslmfrompypr.stattestimport*# Load datarg=genfromtxt('sp.dat')year=rg[:,0].astype(int)sp=rg[:,1]d=11# Number of inputstrain_until_year=1920# includinglast_train=nonzero(year==train_until_year)[0][0]-dX=lm.create_lag_matrix(np.c_[sp],d)# Training and test setsXtrain=X[0:last_train+1,:]Xtest=X[last_train+1:,:]cluster_init_kw={'cluster_init':'kmeans','max_init_iter':5, \
'cov_init':'var','verbose':True}cen_lst,cov_lst,p_k,logL=gmm.em_gm(Xtrain,K=5,max_iter=200, \
delta_stop=1e-5,init_kw=cluster_init_kw,verbose=True)# PredictT=np.c_[X[:,-1]]Y=np.zeros(T.shape)foriinrange(len(Y)):cond_input=np.concatenate((X[i,:-1],np.nan*np.ones(1)))cen_cond,cov_cond,mc_cond=gmm.cond_dist(np.array(cond_input),\
cen_lst,cov_lst,p_k)forjinrange(len(cen_cond)):Y[i]=Y[i]+(cen_cond[j]*mc_cond[j])SE=(Y-T)**2.0# Plot datafigure()subplot(211)plot(year,rg[:,1],'b',label='Target')title('Yearly mean of group sunspot numbers')xlabel('Year')ylabel('Number')plot(year[d:],Y,'r',label='Output')legend()title('Sunspot predcition using a NN with %d inputs'%d)axvline(train_until_year)subplot(212)plot(year[d:],SE)xlabel('Year')ylabel('$(y(x^n)-t^n)^2$')title('Mean square error = %.2f'%np.mean(SE))# Plot sample autocorrelation of the residual of the fitted modelfigure()lags=range(0,21)stem(lags,sac(Y-T,lags))title('Sample Autocorrelation of residuals')xlabel('Lag')ylabel('Sample Autocorrelation')# Predict the future (pass outputs from the ANN to the inputs)curr_input=Xtest[0,:-1]res=[Ttest[0,0]]foriinrange(1,len(Ttest)):printcurr_inputcond_input=np.concatenate((curr_input,np.nan*np.ones(1)))cen_cond,cov_cond,mc_cond=gmm.cond_dist(np.array(cond_input),\
cen_lst,cov_lst,p_k)pred=0forjinrange(len(cen_cond)):pred=pred+(cen_cond[j]*mc_cond[j])curr_input=np.concatenate((curr_input[1:],pred.flatten()))res.append(pred)figure(figsize=(8,3))plot(year[-len(res):],res,'r--',label='Predicted')plot(year,rg[:,1],'b',label='Target')axvline(year[-len(res)])xlabel('Year')ylabel('Number')