what's wrong with my model?

[Fmajor]

Hi:
I'm a beginner in learning MCMC and I only know some concept of Metropolis step methods(I take it as random walk with a multivariate normal distribution)
I notice that I must define ALL the variables as some distribution, for example, I want to fit a function with some parameters but i know nothing about them, so I define them with Uniform('someParameter',lower = MyEstimation, upper = MyEstimation), and define MyModel = Normal('MyModel', mu=functionWithParameters, tau = Observe_data_stds, value = Observe_datas, observed=True), then use M=MCMC('MyModel') and begin to fit, but the trace of the parameters are just like "WWWWWWWW" and the hist "IIIIIII". They are not stable at all and seem to be random sampled from the Uniform distribution they are defined as. Of course the chain can not satisfy the Normal distribution, what's wrong with my model?

-----------------------------------mymodel.py-------------------------------------

#Omiga_m, Omiga_lamuda, H are the three parameters
#SN_z, SN_y, SN_std are (580,) numpy arrays
#SN_z is the independent variable, SN_y the observed data, SN_std the std error of SN_y
Omiga_m = Uniform('Omiga_m',lower = 0, upper = 1)
Omiga_lamuda = Uniform('Omiga_lamuda',lower = 0, upper = 1)
H = Uniform('H',lower = 60, upper = 80)# they must have a prior??????

@deterministic(plot=False)
def miu_th(Omiga_m = Omiga_m, Omiga_lamuda = Omiga_lamuda, H=H):
    return SN_miu_th(SN_z,Omiga_m,Omiga_lamuda,H)
# the function miu_th takes in the three parameter and the independent variable, return a (580,) numpy array
# the likelihood of the parameters should be exp^(-0.5*sum (((miu_th-SN_y)/SN_std)^2))
# so I write the model as 
SN = Normal('SN', mu=miu_th, tau = SN_std, value = SN_y, observed=True)
#miu_th return the 580 y-value calculated, and SN_y is the obsvered data, it has the std error SN_std, which are also (580,) numpy array
#and what should I do if I have the cov matrix of the observed data rather than their own std???

#will, since I can make my model to be 580 Normal distribution,
#can I make the sum((miu_th-SN_y)**2/SN_std**2) to be a chi2 distribution??
#if I can do this, which is faster?

Oh, there is another question about the EXAMPLE StraightLineFit(https://github.com/pymc-devs/pymc/wiki/StraightLineFit)
(I can't get the right result to fit the StraightLineFit by pymc.MCMC and
M.sample(iter=10000,burn=2000))

# this definition
@stochastic
def theta(value=array([2.,5.])):
    """Slope and intercept parameters for a straight line.
    The likelihood corresponds to the prior probability of the parameters."""
    slope, intercept = value
    prob_intercept = uniform_like(intercept, -10, 10)
    prob_slope = log(1./cos(slope)**2)
    return prob_intercept+prob_slope
#is it a custom prior of the parameters with a union distribution??
#so it return the log-probability density?

# is it because the variable x have a std so we define it also as a stochastic??
# should i define my SN_z the same way??
x = Uniform('x', lower=0, upper=50, value=init_x)

#and this
y = Normal('y', mu=modelled_y, tau=2, value=data_y, observed=True)
# that means the observed data have the same std 2?

Is it always right to define the parameters with uniform distributions when I know nothing about it???
will this has a bad effect on the MCMC process??
can I define the cov matrix of the Metropolis step methods(I think it as a random walk with a multivariate normal distribution, so it can have a cov matrix?)

Oh, sorry for SO many questions.....I need to do some work Using the MCMC, and I actually have done it by numpy only and have a not bad result, but now i want to try Pymc, I'm eager to understand it

------------update at 3pm----------
I have tried to change the parameter of pymc to fit my model, but I still can't understand how it works

the posterior distribution of the 3 parameters is "monoclonal peak"
and the union distribution of the 2 parameters Omiga_m, Omiga_lamuda is just like a bivariate normal distribution
before that, I finished my work using numpy like this:

for i in range(length-1):  #begin the MCMC
    print 'i=',i
    
    _lastChain = chain[i,0:3] # it record the last node's 3 parameters Omiga_m, Omiga_lamuda and H
    
    _lastkai = chain[i,3] # this kai(my pool pronunciation....) is the chi2 of 
    # sum((min_th-SN_y)**2/SN_std**2), which contains the information of the posterior's probability density
    
    _newChain = _lastChain + delta_v() # random walk with a multivariate normal distribution
    # the cov matrix of the distribution is the inverse of the Fisher matrix in the best point times a scale
    # I tried the scale many times to make sure the MCMC have a medium acception rate
    
    _newkais = SN_kai(_newChain[0], 1-_newChain[0]-_newChain[1], _newChain[2])(SN_z,SN_y)
    _newkai = np.sum(_newkais**2/SN_std**2)
    # calculate the new chi2
    
    if _newkai>_lastkai:# if the new chi2 is bigger(the probability is smaller 
    # since the posterior's probability density function is like exp(-0.5*chi2))
    # we accept it with a probability of (thisNode's probability density / lastNode's probability density)
    # that is exp(-0.5*(_newkai-_lastkai))
        delta_kai = _newkai-_lastkai
        f = np.exp(-0.5*delta_kai) 
        if f>np.random.rand(): # accept it! make the new node with the new parameters
            chain[i+1,0:3] = _newChain
            chain[i+1,3] = _newkai
            accepted[i+1] = 1
            continue
        else:  # reject it, make the new node witn old parameters
            chain[i+1,0:3] = _lastChain
            chain[i+1,3] = _lastkai
            continue
            
    else: # the new chi2 is smaller (the probability is bigger), accept it at once
        chain[i+1,0:3] = _newChain
        chain[i+1,3] = _newkai
        accepted[i+1] = 1
        continue

this program gives me a not bad result
i have the acception rate of about 40%
the chain of the parameters vibrate near the best point
and the hists are like Normal distributions

while in Pymc
the chain of the parameters vibrate violently(actually the plot range is full with the chain....)
the hists are not Normal-like

I don't understand how Pymc 'tune' the std, can I set it stable just like the above program?
Or I set the wrong args to fail my work?

thank you for your answer~

[fonnesbeck]

It would be useful if you posted your model. Its difficult to diagnose the problem from what you have posted.

[Fmajor]

thank you for your answer :)
I have update my issue, and add some ....more questions

[Fmajor]

I am waiting for your answer @_@...., @fonnesbeck

[fonnesbeck]

You want to let PyMC automatically tune the proposal so that it can achieve a good acceptance rate. Its best not to interfere with this, particularly for a beginner. If you read the user's guide, it gives you some details about this.

Modeling your parameters as a Uniform random variable works in some situations, yes. You have to make sure that it does not exclude regions of the state space that might be part of the support of the true variable.

Your model looks okay, but I dont know what your function SN_miu_th is doing. I suggest running the Python debugger and seeing if you are getting sensible values out of your function at each step. Also, I recommend setting initial values to your stochastic parameters, incase the problem is poorly-chosen starting values.