Does spotpy's Dream not exactly Vrugt, J. A. (2016) proposed Dream?

[BeingHapppy]

Hi @thouska and @bees4ever

Thank you very much for your work. I've been using spotpy for almost three years. Today, I find some differences between spotpy's Dream method and the Dream method proposed by Vrugt, J. A. I'm curious why they are different.

（1）In spotpy, it seems we use only one couple chains to get_new proposal vector, but Vrugt seems use random number of chains. here is the code from spotpy:

` def _get_gamma(self,N):

N = Number of parameters

    p = np.random.uniform(low=0,high=1)
    if p >=0.2:
        gamma = 2.38/np.sqrt(2*int(N))#/self.gammalevel
    else:
        gamma = 1
    return gamma

def get_other_random_chains(self,cur_chain):
    valid=False        
    while valid == False:         
        random_chain1 = np.random.randint(0,self.nChains)
        random_chain2 = np.random.randint(0,self.nChains)
        if random_chain1!=cur_chain and random_chain2!=cur_chain and random_chain1!=random_chain2:
            valid=True
    return random_chain1, random_chain2

def get_new_proposal_vector(self,cur_chain,newN,nrN):
    gamma = self._get_gamma(nrN)
    random_chain1,random_chain2 = self.get_other_random_chains(cur_chain)
    new_parameterset=[]        
    #position = self.chain_samples-1#self.nChains*self.chain_samples+self.chain_samples+cur_chain-1
    cur_par_set = list(self.bestpar[cur_chain][self.nChainruns[cur_chain]-1])
    random_par_set1 = list(self.bestpar[random_chain1][self.nChainruns[random_chain1]-1])
    random_par_set2 = list(self.bestpar[random_chain2][self.nChainruns[random_chain2]-1])

    for i in range(self.N):#Go through parameters

        if newN[i] == True:
            new_parameterset.append(cur_par_set[i] + gamma*np.array(random_par_set1[i]-random_par_set2[i]) + np.random.normal(0,self.eps))
        else:
            new_parameterset.append(cur_par_set[i])

    new_parameter=self.check_par_validity_reflect(new_parameterset)        
    #new_parameter=self.check_par_validity_bound(new_parameterset)        
    return new_parameter`

below is from Vrugt:

The highlighted part in the figure is the evidence that Vrugt uses multiple chains.

And Vrugt also point out the total number of chains should equal 2*sita+1, sita is the max number of couple chains when proposal vetor.

(2) outlier chain detection

In the paper of Vrugt2016, he says that he use outlier chain detection during the sampling, but he did not show how to do. It seems spotpy Dream does not have outlier chain detection.

[thouska]

Hi @BeingHapppy, thank you for your detailed message and also for using spotpy for such a long time! Your describtion looks very valid and is as such a feature that is missing in spotpy, so far. May I ask: Do you think you could provide us a pull request, which implents the missind features in Dream? I can help on the way, if needed.

[BeingHapppy]

Hi @thouska, thank you for your reply! Implementing the missing features is a challenge for me, but I'll try. I appreciate that you said you can provide help.

[bees4ever]

Thanks @BeingHapppy for reporting. Unfortunately as @thouska also mentioned, for concrete implementation I also currently don't have time. Let's see what's going on in autumn. Cool that you use this tool.

[BeingHapppy]

Thank you for your reply. MCMC method is the foundation of our team's research, and spotpy dream is the easiest tool we've ever used. Before meeting spotpy, we spent almost a whole year learning MCMC and still couldn't program it. I just did very little work and thank you for opening the doors of MCMC to us. Looking forward to the future of MCMC, your work is amazing.

Best wishes, Yantong

------------------ Original message ------------------ From: "Ben"; Sendtime: Wednesday, Jun 22, 2022 4:31 AM To: "thouska/spotpy"; Cc: "Yantong @.***>; "Mention"; Subject: Re: [thouska/spotpy] Does spotpy's Dream not exactly Vrugt, J. A. (2016) proposed Dream? (Issue #284)

Thanks @BeingHapppy for reporting. Unfortunately as @thouska also mentioned, for concrete implementation I also currently don't have time. Let's see what's going on in autumn. Cool that you use this tool.

— Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: @.***>