ws$accept hits 0 after first iteration again (for parallelization)

[gtng92]

In commit 926609f, ran a simulation using parallelization implementation according to issue #26. After about 20 minutes, the following error:

run.smc niter:  1 
ws$epsilon:  0.8270003 
config$final.epsilon:  0.1 
result$accept.rate:  0 
config$final.accept.rate:  0.1 
Error in density.default(temp$beta, weights = temp$weight) : 
  need at least 2 points to select a bandwidth automatically
Calls: lines -> density -> density.default
In addition: Warning message:
In density.default(temp$beta, weights = temp$weight) :
  sum(weights) != 1  -- will not get true density

The error is is due to plotting, but the main concern is hitting the accept rate at the first iteration. Looking into the mh.ratio's again.

[gtng92]

Sidenote: also tried this with Mathias's Yule model and result$accept.rate also hits 0 after first iteration.

[gtng92]

Switched to checking if parallelization works for a simpler model (yule).

old.nbhd <- sum(ws$dists[,i] < ws$epsilon)    # how many samples are in the neighbourhood of data?
new.nbhd <- sum(new.dists[,i] < ws$epsilon)
mh.ratio <- mh.ratio * new.nbhd / old.nbhd

Variables old.nbhd and new.nbhd sometimes output value of 0, which turns the mh.ratio into NaN. This is stochastic. This did not happen before the parallelization implementation.

Otherwise parallelization seems to work in shortening runtime.

[gtng92]

I added two clauses to check for NaN, as well as if there is a returned NULL (from rejected particles)

if (length(i) == 0) {     
    next
}
if (NaN %in% i) {
    next
}

Hopefully this is fine?

[ArtPoon]

We'll have to test two versions (with and without these next clauses) on simulated data to see how this affects prediction accuracy.

[ArtPoon]

I think the problem is how the new.dists object is being returned from mclapply. Pass the row vector new.dists[i] into the return objects that are gathered as the return value (list) for mclapply, and then reconstitute a matrix ordered by row indices (first value in return object) outside of mclapply.

[gtng92]

I passed new.dists as a vector in the return object rather than populating a matrix of new.dists outside of the mclapply, and I compared parallelization to check that merging results from threads is working properly. I copied a sample:

n	part.num	weight	lambda	dist.1	dist.2	dist.3	dist.4	dist.5
1	1	0.0118483412	0.30315	0.15673	0.20844	0.1943	0.21388	0.18314
1	2	0.009478673	1.09258	0.29435	0.29231	0.32897	0.29269	0.29372
1	3	0.009478673	1.73045	0.30696	0.33262	0.30416	0.32755	0.33092
1	4	0.009478673	1.12879	0.33567	0.29004	0.31124	0.29125	0.30322
1	5	0.0118483412	2.68067	0.31435	0.32743	0.33046	0.31893	0.31246
1	6	0.0118483412	0.66457	0.27805	0.27929	0.26469	0.27106	0.27349
1	7	0.0118483412	1.22272	0.29534	0.29042	0.30667	0.30001	0.30571
1	8	0.0118483412	0.19587	0.16099	0.1386	0.13418	0.10999	0.15808
1	9	0.0071090047	2.12755	0.3267	0.34033	0.34465	0.31512	0.31059
1	10	0.0118483412	1.19267	0.29325	0.29585	0.32507	0.29361	0.2883
1	11	0.0118483412	1.36457	0.30178	0.31357	0.29698	0.30055	0.30156
1	12	0.0118483412	1.27433	0.29998	0.32996	0.31197	0.31464	0.3025

When perturbing particles, if the particle was accepted, the particle was written to a separate table. Below particle 10 and 11 were rejected, and the rest were retained and incorporated into ws$dists

part.num	mh.ratio	lambda	dist.1	dist.2	dist.3	dist.4	dist.5
1	1.0043417079	0.3031458446	0.1567258973	0.2084382774	0.1942950921	0.2138820201	0.183136695
2	1.255850512	1.0925795619	0.2943499709	0.2923068338	0.3289730025	0.2926916396	0.2937162621
3	1.2050725618	1.7304483532	0.3069638106	0.3326150027	0.3041640263	0.3275458928	0.3309234589
4	0.996839063	1.1287941702	0.3356732582	0.2900407844	0.311241625	0.2912462919	0.3032218013
5	1.0037199224	2.6806748366	0.3143498225	0.3274334579	0.3304630712	0.3189295192	0.3124646548
6	1.0948451332	0.6645661761	0.2780547839	0.2792937201	0.264685498	0.2710595922	0.2734917836
7	1.0124995482	1.222719962	0.2953396574	0.2904201242	0.3066665508	0.3000068714	0.305710307
8	2.2772833491	0.1958671613	0.1609945419	0.1386034062	0.1341826656	0.109986818	0.1580814178
9	0.9683790581	2.1275542956	0.3267005721	0.3403317498	0.3446453467	0.3151235459	0.3105892557
12	0.989995566	1.2743262198	0.299983392	0.3299570182	0.3119710161	0.3146444966	0.302502847

[gtng92]

I also removed the if clause dealing with the NaN error, and running smc worked for 100 particles. Running Yule model with 1000 particles to double check, and closing after.