建模杂谈系列253 序列突变点的判定

2024-09-03 01:04

本文主要是介绍建模杂谈系列253 序列突变点的判定,希望对大家解决编程问题提供一定的参考价值,需要的开发者们随着小编来一起学习吧!

说明

使用pycm3进行推断。

内容

1 环境搭建

使用conda创建对应的包环境,然后再通过jupyter运行

conda create -c conda-forge -n pymc_env "pymc>=5"
conda activate pymc_envpip3 install  ipython -i https://mirrors.cloud.tencent.com/pypi/simple
pip3 install  jupyter -i https://mirrors.cloud.tencent.com/pypi/simplenohup jupyter lab --allow-root --ip='*' --port=8888 > /dev/null 2>&1 &

嗯,我发现启动jupyter后竟然没有pymc,又重装了一下。最近阿里的镜像慢的不行了,腾讯快。
!pip3 install pymc -i https://mirrors.cloud.tencent.com/pypi/simple

2 实验

判定一个时间序列,产生突变点的位置

主要算是重启一下之前的记忆,确保新的环境搭建成功

%matplotlib inline
from IPython.core.pylabtools import figsize
import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl
mpl.style.use("ggplot")
# figsize(11, 9)figsize(12.5, 3.5)
count_data = np.loadtxt("txtdata.csv")
n_count_data = len(count_data)
plt.bar(np.arange(n_count_data), count_data, color="#348ABD")
plt.xlabel("Time (days)")
plt.ylabel("count of text-msgs received")
plt.title("Did the user's texting habits change over time?")
plt.xlim(0, n_count_data);

说的是作者自己发短信的数据,有一天发生了变化,但是从图上看不出来。
在这里插入图片描述

然后构建了一个模型框架,假设有两种不同的模式

import pymc as pmwith pm.Model() as model:alpha = 1.0/count_data.mean()  # Recall count_data is the# variable that holds our txt countslambda_1 = pm.Exponential("lambda_1", alpha)lambda_2 = pm.Exponential("lambda_2", alpha)tau = pm.DiscreteUniform("tau", lower=0, upper=n_count_data - 1)with model:idx = np.arange(n_count_data) # Indexlambda_ = pm.math.switch(tau > idx, lambda_1, lambda_2)with model:observation = pm.Poisson("obs", lambda_, observed=count_data)

然后将观察数据给到模型,执行mcmc,然后画图

### Mysterious code to be explained in Chapter 3.
with model:step = pm.Metropolis()trace = pm.sample(10000, tune=5000, step=step, return_inferencedata=False)lambda_1_samples = trace['lambda_1']
lambda_2_samples = trace['lambda_2']
tau_samples = trace['tau']figsize(12.5, 10)
#histogram of the samples:ax = plt.subplot(311)
ax.set_autoscaley_on(False)plt.hist(lambda_1_samples, histtype='stepfilled', bins=30, alpha=0.85,label="posterior of $\lambda_1$", color="#A60628", density=True)
plt.legend(loc="upper left")
plt.title(r"""Posterior distributions of the variables$\lambda_1,\;\lambda_2,\;\tau$""")
plt.xlim([15, 30])
plt.xlabel("$\lambda_1$ value")ax = plt.subplot(312)
ax.set_autoscaley_on(False)
plt.hist(lambda_2_samples, histtype='stepfilled', bins=30, alpha=0.85,label="posterior of $\lambda_2$", color="#7A68A6", density=True)
plt.legend(loc="upper left")
plt.xlim([15, 30])
plt.xlabel("$\lambda_2$ value")plt.subplot(313)
w = 1.0 / tau_samples.shape[0] * np.ones_like(tau_samples)
plt.hist(tau_samples, bins=n_count_data, alpha=1,label=r"posterior of $\tau$",color="#467821", weights=w, rwidth=2.)
plt.xticks(np.arange(n_count_data))plt.legend(loc="upper left")
plt.ylim([0, .75])
plt.xlim([35, len(count_data)-20])
plt.xlabel(r"$\tau$ (in days)")
plt.ylabel("probability");

最后可以推断在45天的时候发生跳变(作者自己承认了这一点)
在这里插入图片描述

3 总结

  • 1 这种模型可以推断看不见的事实
  • 2 matplotlib用好了也可以画出很好看的图

注意:理论上,pymc是可以借用显卡计算的。只不过之前它的后端比较奇葩,现在可能失传了。然后我不太清楚的是现在怎么和显卡挂上,单靠cpu只能做一些小规模的计算。

这篇关于建模杂谈系列253 序列突变点的判定的文章就介绍到这儿,希望我们推荐的文章对编程师们有所帮助!



http://www.chinasem.cn/article/1131577

相关文章

Spring Security 从入门到进阶系列教程

Spring Security 入门系列 《保护 Web 应用的安全》 《Spring-Security-入门(一):登录与退出》 《Spring-Security-入门(二):基于数据库验证》 《Spring-Security-入门(三):密码加密》 《Spring-Security-入门(四):自定义-Filter》 《Spring-Security-入门(五):在 Sprin

科研绘图系列:R语言扩展物种堆积图(Extended Stacked Barplot)

介绍 R语言的扩展物种堆积图是一种数据可视化工具,它不仅展示了物种的堆积结果,还整合了不同样本分组之间的差异性分析结果。这种图形表示方法能够直观地比较不同物种在各个分组中的显著性差异,为研究者提供了一种有效的数据解读方式。 加载R包 knitr::opts_chunk$set(warning = F, message = F)library(tidyverse)library(phyl

【生成模型系列(初级)】嵌入(Embedding)方程——自然语言处理的数学灵魂【通俗理解】

【通俗理解】嵌入(Embedding)方程——自然语言处理的数学灵魂 关键词提炼 #嵌入方程 #自然语言处理 #词向量 #机器学习 #神经网络 #向量空间模型 #Siri #Google翻译 #AlexNet 第一节:嵌入方程的类比与核心概念【尽可能通俗】 嵌入方程可以被看作是自然语言处理中的“翻译机”,它将文本中的单词或短语转换成计算机能够理解的数学形式,即向量。 正如翻译机将一种语言

uva 10131 最长子序列

题意: 给大象的体重和智商,求体重按从大到小,智商从高到低的最长子序列,并输出路径。 代码: #include <iostream>#include <cstdio>#include <cstdlib>#include <algorithm>#include <cstring>#include <cmath>#include <stack>#include <vect

poj 1127 线段相交的判定

题意: 有n根木棍,每根的端点坐标分别是 px, py, qx, qy。 判断每对木棍是否相连,当他们之间有公共点时,就认为他们相连。 并且通过相连的木棍相连的木棍也是相连的。 解析: 线段相交的判定。 首先,模板中的线段相交是不判端点的,所以要加一个端点在直线上的判定; 然后,端点在直线上的判定这个函数是不判定两个端点是同一个端点的情况的,所以要加是否端点相等的判断。 最后

基于UE5和ROS2的激光雷达+深度RGBD相机小车的仿真指南(五):Blender锥桶建模

前言 本系列教程旨在使用UE5配置一个具备激光雷达+深度摄像机的仿真小车,并使用通过跨平台的方式进行ROS2和UE5仿真的通讯,达到小车自主导航的目的。本教程默认有ROS2导航及其gazebo仿真相关方面基础,Nav2相关的学习教程可以参考本人的其他博客Nav2代价地图实现和原理–Nav2源码解读之CostMap2D(上)-CSDN博客往期教程: 第一期:基于UE5和ROS2的激光雷达+深度RG

flume系列之:查看flume系统日志、查看统计flume日志类型、查看flume日志

遍历指定目录下多个文件查找指定内容 服务器系统日志会记录flume相关日志 cat /var/log/messages |grep -i oom 查找系统日志中关于flume的指定日志 import osdef search_string_in_files(directory, search_string):count = 0

POJ1631最长单调递增子序列

最长单调递增子序列 import java.io.BufferedReader;import java.io.InputStream;import java.io.InputStreamReader;import java.io.PrintWriter;import java.math.BigInteger;import java.util.StringTokenizer;publ

数学建模笔记—— 非线性规划

数学建模笔记—— 非线性规划 非线性规划1. 模型原理1.1 非线性规划的标准型1.2 非线性规划求解的Matlab函数 2. 典型例题3. matlab代码求解3.1 例1 一个简单示例3.2 例2 选址问题1. 第一问 线性规划2. 第二问 非线性规划 非线性规划 非线性规划是一种求解目标函数或约束条件中有一个或几个非线性函数的最优化问题的方法。运筹学的一个重要分支。2

leetcode105 从前序与中序遍历序列构造二叉树

根据一棵树的前序遍历与中序遍历构造二叉树。 注意: 你可以假设树中没有重复的元素。 例如,给出 前序遍历 preorder = [3,9,20,15,7]中序遍历 inorder = [9,3,15,20,7] 返回如下的二叉树: 3/ \9 20/ \15 7   class Solution {public TreeNode buildTree(int[] pr