bnlearn::bn.fit difference and calculation of methods "mle" and "bayes"
I try to understand the differences between the two methods bayes
and mle
in the bn.fit
function of the package bnlearn
.
I know about the debate between the frequentist and the bayesian approach on understanding probabilities. On a theoretical level I suppose the maximum likelihood estimate mle
is a simple frequentist approach setting the relative frequencies as the probability. But what calculations are done to get the bayes
estimate? I already checked out the bnlearn documenation, the description of the bn.fit function and some application examples, but nowhere there's a real description of what's happening.
I also tried to understand the function in R by first checking out bnlearn::bn.fit
, leading to bnlearn:::bn.fit.backend
, leading to bnlearn:::smartSapply
but then I got stuck.
Some help would be really appreciated as I use the package for academic work and therefore I should be able to explain what happens.
See also questions close to this topic

Adjust the height of NavBar menu in a shiny app
I have a basic shiny app in which I want to be able to set the height of the Navbar menu named Navnar! I do not see any choice for setting its height.
#ui.r library(markdown) navbarPage("Navbar!", tabPanel("Plot", sidebarLayout( sidebarPanel( ), mainPanel( ) ) ), tabPanel("Summary" ), navbarMenu("More", tabPanel("Table" ) ) ) #server.r function(input, output, session) { }

How to use monte carlo to plot histogram?
enter image description hereI am trying to run a monte carlo simulation but don't know how to write the code. I need to plot a distribution of rolling sixes in 10 rolls of a die without using rbinom.
I have tried using the replicate function but my histogram does not look right.
dat < replicate(n=10, expr = sum(runif(10, min=0, max=1)>1/6)) hist(dat, prob=TRUE)
I think since it is a simulation I am expected to have 10,000 trials (according my professor) but I don't know where to put the number of trials.enter image description here

Language R: Variable doesn't change inside a 'for()' cycle?
I am trying to analyze the value of "edoActual", but it doesn't seem to change its value after the for cycle.
When running it in RGui, it doesn't point any errors or warnings. When running it on RStudio it goes nuts. Either way it doesn't work.
edoActual = "q0"; print("Input cadena: "); cadena < scan("", what="char"); for(i in 1:length(cadena)) { if( edoActual == "q0" && cadena[i] == "a" ) { edoActual = "q1"; } if( edoActual == "q0" && cadena[i] == "b" ) { edoActual = "q0"; } if( edoActual == "q1" && cadena[i] == "a" ) { edoActual = "q0"; } if( edoActual == "q1" && cadena[i] == "b" ) { edoActual = "q1"; } if( cadena[i] != "a" && cadena[i] != "b") { print("Error: invalid."); break; } } if( i == length(cadena) ) { if(edoActual == "q1") {print("accepted (:");} if(edoActual == "q0") {print("rejected ):");} }
What is supposed to happen:
[1] "Input cadena: " 1: a 2: a 3: b 4: a 5: Read 4 ítems [1] "accepted (:"

Can anyone help me with adding a slump into rbacon? It looks simple but it keeps throwing an error
I'm using rbacon to create an age model. I am using the default core that comes with the package "MSB2K". The documentation says to add a slump, insert the code: slump=c(). Requires pairs of depths, e.g., slump=c(10,15,60,67) for slumps at 6760 and 1510 cm core depth.
I try the following code but it gives me an error
Bacon("MSB2K", slump = c(30, 34))
Error in if (!is.na(hiatus.depths)[1]) { : missing value where TRUE/FALSE needed
Does anyone know why this might be?

How to plot trace of the logposterior (up to proportionality)?
I am working to get some trace plots of posterior from a JAGS model. What is the best way to get the trace plot of the logposterior (up to proportionality)?
mydata = read.table('aaa', header = F) jmodel = jags.model(file = '111.model', data = list(x = mydata, N = 50, k=2), n.chains = 2) samps_f = coda.samples(jmodel, variable.names=c('a'), n.iter=1e3)

building PyMC3 model incorporating different measurements
I am trying to incorporate different types and replicates of measurements into one model in PyMC3.
Consider the following model: P(t)=P0*exp(kBt) where P(t), P0, and B are concentrations. k is a rate. We measure P(t) at different times and B once, all through counting of particles. k is the parameter of interest we are trying to infer.
My question has two parts: (1) How to incorporate measurements on P(t) and B into one model? (2) How to use a variable number of replicate experiments to inform on the value of k?
I think I can answer part (1), but am unsure about whether it is right or done in the right flavour. I failed to generalise the code to include a variable number of replicates.
For one experiment (one replicate):
times=np.asarray([time0,time1,...]) counts=np.asarray([countforB,countforPattime0,countforPattime1,...]) basic_model = pm.Model() with basic_model: k=pm.Uniform('k',0,20) B=pm.Uniform('B',0,1000) P=pm.Uniform('P',0,1000) exprate=pm.Deterministic('exprate',k*B) modelmu=pm.math.concatenate(B*(np.asarray([1.0]),P*pm.math.exp(exprate*ts))) Y_obs=pm.Poisson('Y_obs',mu=modelmu,observed=counts))
I tried to include different replicates along the lines of the above, but to no avail:
... k=pm.Uniform('k',0,20) # same within replicates B=pm.Uniform('B',0,1000,shape=numrepl) # can vary between expts. P=pm.Uniform('P',0,1000,shape=numrepl) # can vary between expts. exprate=??? modelmu=???

Integrate weight to a Bayesian network modal
I want create a Bayesian network. But my problem is, I can't integrate the weighting of my data in the Bayesian network (For my Bayesian Network I use package bnlearn, available with R).
Please somebody can help me to find a solution!
For exemple:
test = head(colon[,c(4,10)]) test$poids = c(0.2,0.6,0.5,0.6,0.9,0.8)
Result 1
Without weighting, I have :
table(test$sex,test$status)
Result 2
And when I use weight, Ihave:
wtd.table(test$sex,test$status, weight = test$poids)
I want the same result (result 2), but using Bayesian network.

How to perfrom MAP estimation in a Bayesian Network in R?
I have to perform the Bayesian Network test, using MAP estimation (Maximum A Posteriori). I have a (gaussian) dataframe consisting of 5 variables, varying from 2 till 16 unique values per variable. Does anyone know who to set the arcs using the bnlearn package and how to perform the estimate learning? Right now I receive the following error:
Error in check.fitting.method(method, data) : Bayesian parameter estimation for (conditional) Gaussian Bayesian networks is not implemented.
#MAP estimate learning library(bnlearn) pdag = iamb(BayesianDF) pdag #How to set arcs for the other nodes? dag = set.arc(pdag, from = "results.cluster", to = "type_touch") fit = bn.fit(dag, BayesianDF, method = "bayes") fit
I expect to get parameter estimates which represent the (conditional) probabilities of the Bayesian Network.

Prediction model for bayesian network in Rjags
I have a text to predict loglevel of a poisonous gas in houses in different counties in Minnesota. This is the part of the table I have been given:
The whole table is located here:
Now, the loglevel of uranium depends on loglevel of uranium in the soil (uu), whether the house has a basement (floor=1 is with basement, floor=0 is without) and in what county the house is based. Hence I came up with this hierarchical model:
Now, the task is "Predict the logradon level of a house with a basement for each of the 84 counties" in the dataset. And I have to use rjags to do it. The only problem is, I have no idea how to implement this model. Could someone please lend me a hand? I took a course in bayesian statistics, and I am really really bad with R. On top of that the material is quite hard. Any help will be much appreciated.

Running maximum likelihood estimation using mle2 returns error: subscript out of bounds
I am trying to do a maximum likelihood estimation of some parameter values for an ordinary differential equation (ODE) model. I have used optim() without a problem, but I would like to instead use "mle2" from the package bbmle (because of ease for calculating confidence intervals).
While the optim() command runs fine, mle2() returns an error
Error in solution.prev[2:5, "op_iggP"] : subscript out of bounds
Apologies for the length of the code posted. For it to be reproducible, I have to include this much.
Any help much appreciated.
# Reproducible example ## ======== Main model ======== ## library(deSolve) ab_mod = function(age, state, parms) { with(as.list(c(age, state, parms)), { ## Timevarying lambda function l_exp_damp = lambda0 * (age * exp(gradient*age)) foi < l_exp_damp # ODEs dtp_iggP < ( foi * tp_iggN )  iggR * tp_iggP dtp_iggN < ( foi * tp_iggN ) + iggR * tp_iggP dtp_tgerpP< (p1* foi) * tp_tgerpN  tgerpR * tp_tgerpP dtp_tgerpN< (p1*foi) * tp_tgerpN + tgerpR * tp_tgerpP op_iggP < tp_iggP * (se_igg + sp_igg 1) + (1sp_igg ) op_tgerpP_roc < tp_tgerpP * (se_tgerp_roc + sp_tgerp_roc 1) + (1sp_tgerp_roc ) op_tgerpP_sd < tp_tgerpP * (se_tgerp_sd + sp_tgerp_sd 1) + (1sp_tgerp_sd ) return(list(c( dtp_iggP, dtp_iggN, dtp_tgerpP, dtp_tgerpN), op_iggP=op_iggP, op_tgerpP_roc=op_tgerpP_roc, op_tgerpP_sd=op_tgerpP_sd, foi=foi )) }) } # === Define parameter values for solving model ==== # age = c(0,3.5,11,18,25) state = c(tp_iggP=0, tp_iggN=1, tp_tgerpP=0, tp_tgerpN=1) P < c(lambda = 0.08, iggR = 0, tgerpR = 2, se_igg = 0.95, sp_igg = 0.92, se_tgerp_roc=0.915, sp_tgerp_roc = 0.754, se_tgerp_sd= 0.85, sp_tgerp_sd = 0.95, p1 = 1.0, gradient = 0.05, lambda0 = 0.02) # =================== Solve model =================== # solution = ode(y = state, times = age, parms = parms, func = ab_mod) # =================== Create data =================== # data1<matrix(nrow=4, ncol=5) data1[1,] < c(3.5, 13, 9, 36, 116) data1[2,] < c(11, 59, 20, 75, 156) data1[3,] < c(18, 33, 13, 34, 48) data1[4,] < c(25, 45, 31, 45, 60) colnames(data1) < c("age","k_igg", "k_tgerp_sd", "k_tgerp_roc", "n") data1<as.data.frame(data1) # ============ Joint likelihood function ============ # loglik < function(kIgG, nIgG, prevIgG, kTgERProc, nTgERProc, prevTgERProc, kTgERPsd, nTgERPsd, prevTgERPsd){ dbinom(kIgG, nIgG, prevIgG, log=T) + dbinom(kTgERProc, nTgERProc, prevTgERProc, log=T) } # =========== Likelihood wrapper function =========== # loglik.wrap < function(par, data){ # Parameters that I want to estimate P["lambda0"] < par[1] P["gradient"] < par[2] P["tgerpR"]<par[3] # Call ODE model solution.prev < ode(y = state, times = age, parms = P, func = ab_mod) # Select data that I want from model output prevIgG < solution.prev[2:5,"op_iggP"] prevTgERProc < solution.prev[2:5,"op_tgerpP_roc"] prevTgERPsd < solution.prev[2:5,"op_tgerpP_sd"] logliks < loglik(kIgG = data1$k_igg, kTgERProc = data1$k_tgerp_roc, kTgERPsd = data1$k_tgerp_sd, nIgG = data1$n, nTgERProc = data1$n, nTgERPsd = data1$n, prevIgG = prevIgG, prevTgERProc = prevTgERProc, prevTgERPsd = prevTgERPsd) return(sum(logliks)) } ## Check that log likelihood function works par < c(0.02, 0.05, 2) mod1<loglik.wrap(par=par, data=data1) mod1 # ==================== MLE analysis ===================== # library(bbmle)
optim method: runs fine
res < optim(par=c(0.02, 0.05, 2), fn=loglik.wrap, data=data1, method="LBFGSB", lower=c(0, 0, 0), upper=c(Inf, 1, 1)) res
mle2 method: returns error "Error in solution.prev[2:5, "op_iggP"] : subscript out of bounds"
par < c(0.02, 0.05, 0.12) parnames<c("lambda0","gradient","tgerpR") (mod < mle2(loglik.wrap, start=list(par), parnames = parnames, method="LBFGSB", data=data1))

Total newbie not into programming but need to learn how to get to this kind of model
(Just a quick reminder that I really have no idea about python and I don't know if person with no experience can do this stuff)
So I'm trying to make a logistic regression model using MLE method to predict the outcome of certain games, I have all the CSV files but I totally can't get how to make this model: Screenshot of model
I get all kind of errors, "NameError: name 'MapWon' is not defined", "TypeError: cannot perform reduce with flexible type".
I literally searched through whole google trying to see how to make those kind of models and there is a lot of things that I don't understand or they don't use exactly MLE method.
Beginning of this code works and it looks like this:
import pandas as pd import numpy as np from sklearn import preprocessing import matplotlib.pyplot as plt plt.rc("font", size=14) from sklearn.linear_model import LogisticRegression from sklearn.model_selection import train_test_split import seaborn as sns sns.set(style="white") sns.set(style="whitegrid", color_codes=True) AGO_data= pd.read_csv("AGO_All.csv") AGO_data= AGO_data.dropna() print(AGO_data.shape) print(list(AGO_data.columns)) AGO_data.head(120)
And it doesn't cause any trouble, but then I have no idea if I have to "train" the data/model or do whatever or can I simply try to make a model that I showed with just csv file?
This is the code I try for strictly this model (I know it's not MLE method but couldn't find anything else):
import statsmodels.api as sm logit_model=sm.Logit(MapWon) result=logit_model.fit() print(result.summary2())
And just getting errors, mostly because "MapWon is not defined"(but it doesn't have problems with commands like:
AGO_data['MapWon'].value_counts()
).And if I use (y,X) just like in tutorial instead of (MapWon) I get "TypeError: cannot perform reduce with flexible type". So is it simply to tough for someone new or if not then what do I do to achieve this model? 
How to plot a weibull3 distribution having all three parameters?
I am fitting my data (wind data) to different distributions (two and three parameter weibull and gamma distributions). I used the method explained here to estimate the parameters of the threeparameter weibull distribution: Weibull Distribution parameter estimation error I used fitdist (from fitdistrplus library) to estimate the parameters of the twoweibull and gamma distribution. Now I want to compare the fit of the twoparameter weibull, threeparameter weibull and gamma distribution to my data on the same plot. I was able to compare the fit of twoparameter weibull and gamma only but failed to include threeparameter weibull. I would highly appreciate any suggestions. Thank you.
estimating parameters of threeparameter and twoparameter weibull:
library(bbmle) library(FAdist) y = my.data eps < 1e8 weibull3 < mle2(y ~ dweibull3(shape=shape, scale = scale, thres = thres), start=list(shape=1, scale=1, thres=5),lower=c(shape=0,scale=0,thres=Inf), upper=c(shape=Inf,scale=Inf,thres=min(y)eps), method="LBFGSB", data=data.frame(y)) weibull3 call: mle2(minuslogl = y ~ dweibull3(shape = shape, scale = scale, thres = thres), start = list(shape = 1, scale = 1, thres = 5), method = "LBFGSB", data = data.frame(y), lower = c(shape = 0, scale = 0, thres = Inf), upper = c(shape = Inf, scale = Inf, thres = min(y)  eps)) coefficients: shape scale thre 2.4828383 5.5248459 0.2530033 Loglikelihood: 1358214 library (fitdistrplus) weibull2 = fitdist(my.data, "weibull") weibull2 Fitting of the distribution ' weibull ' by maximum likelihood Parameters: estimate Std. Error shape 2.324290 0.002344834 scale 5.238753 0.002980665

BNlearn, gRain pEvidence issue with a Bayesian Network
I'm trying to get the likelihood of a data record to belong to a data set characterized by a Bayesian Network.
Here is my current simple code:
library(bnlearn) library(gRain) dag = hc(learning.test) fitted.bnlearn = bn.fit(dag, learning.test, method = "bayes", iss=1) fitted.grain = as.grain(fitted.bnlearn) retractEvidence(fitted.grain) evidence_probability < setEvidence(fitted.grain, c("A","B", "C", "D", "E", "F"), c("Z", "Z", "Z", "Z", "Z", "Z")) evidence_probability pEvidence(evidence_probability)
this outputs:
> pEvidence(evidence_probability) [1] 1
What definitely is not what I'm expecting. What I'm expecting is that by using the smoothing (iss=1) it should output a very small but not null probability.
There is definitely something that I don't understand (all 6 variables have on ly a, b or c values int the "learning.test" training data set).

Probability of a record to belong to the same dataset given a Bayesian Network built with BNlearn
I'm trying to determine the probability of a new record to belong to an existing dataset. I'm using the BNlearn R package to build a Bayesian Network using a large training set.
I then want to assess how anomalous a new record is. For this I want to get a probability for a record for which I have full evidence but don't need to predict any variable.
The pcquery method seems to require at least one variable to predict. The documentation states that the predict method will ignore entries with full evidence.
I spent a day searching the BNlearn documentation without success. So I think it is either not possible with BNlearn or I'm missing the right vocabulary to find what I need in the docs.
Any insights from someone who has more experience with BNlearn is welcome.