karlinMonteCarlo_double: Monte Carlo - Karlin for real scores[p-value]

In localScore: Package for Sequence Analysis by Local Score

Monte Carlo - Karlin for real scores[p-value]

Description

Estimates p-value, for integer scores, based on a Monte Carlo estimation of Gumble parameters from simulations of smaller sequences with same distribution. Appropriate for great sequences with length > 10^3, for i.i.d and markovian sequence models.

Usage

karlinMonteCarlo_double(
  local_score,
  sequence_length,
  simulated_sequence_length,
  FUN,
  ...,
  numSim = 1000,
  plot = TRUE
)

Arguments

local_score	local score observed in a segment.
sequence_length	length of the sequence
simulated_sequence_length	length of simulated sequences produced by FUN
FUN	function to simulate similar sequences with.
...	parameters for FUN
numSim	number of sequences to create for estimation
plot	boolean value if to display plots for cumulated function and density

Details

The length of the simulated sequences is an argument specific to the function provided for simulation. Thus, it has to be provided also in the parameter simulated_sequence_length in the arguments of the "Monte Carlo - Karlin" function. It is a crucial detail as it influences precision and computation time of the result. Note that to get an appropriate estimation, the given average score must be non-positive.

Value

Floating value corresponding to the probability to obtain a local score with value greater or equal the parameter

Examples

new = sample(-7:6, replace = TRUE, size = 1000) 
#MonteCarlo taking random sample from the input sequence itself

karlinMonteCarlo_double(local_score = 66, sequence_length = 1000,  
               FUN = function(x, simulated_sequence_length) {return(sample(x = x, 
               size = simulated_sequence_length, replace = TRUE))}, 
               x=new, simulated_sequence_length = 1000,  numSim = 1000)
      
 #Markovian example (longer computation)
 MyTransMat_reels <-  matrix(c(0.3,0.1,0.1,0.1,0.4, 0.2,0.2,0.1,0.2,0.3, 0.3,0.4,0.1,0.1,0.1, 
0.3,0.3,0.1,0.2,0.1, 0.2,0.1,0.2,0.4,0.1),ncol = 5, byrow=TRUE)

karlinMonteCarlo(local_score = 10.5,sequence_length=200,FUN = transmatrix2sequence, 
matrix = MyTransMat_reels, score =c(-1,-0.5,0,0.5,1),length=1500, 
plot=FALSE, numSim = 1500, simulated_sequence_length =1500)

localScore documentation built on Nov. 3, 2023, 1:08 a.m.