text:
Followers 0

# generic and overloads.

## 4 posts in this topic

having a bit of a problem with something i'd like to do in C#.  in c++ with templates it'd work, since the class would be compiled as w/e type it's being used for.  but in C# it seems that doesn't work in a similar manner.

here's what i'm trying to do: [url]http://ideone.com/D6Tgke[/url]

using System;
namespace Main{

class Test<T>{

public Test(T Value){
X.Foo(Value);
}
}

class X{
public static void Foo(int value){
Console.WriteLine("Foo Int: "+value);
}

public static void Foo(float value){
Console.WriteLine("Foo Float: "+value);
}

public static void Main(){
Test<float> f = new Test<float>(5.0f);
Test<int> i = new Test<int>(5);
return;
}
}
}

which gives me the following errors:

prog.cs(7,15): error CS1502: The best overloaded method match for Main.X.Foo(int)' has some invalid arguments
prog.cs(13,28): (Location of the symbol related to previous error)
prog.cs(7,15): error CS1503: Argument #1' cannot convert T' expression to type int'
Compilation failed: 2 error(s), 0 warnings

any tips on what i could do to achieve this?

0

##### Share on other sites
C# generics are really cumbersome to use with value types.  As far as I know you have to move the type resolution from compile time to runtime. You could make an overload of X.Foo which takes an object and then calls the appropriate method using a series of "Value is Type" if/else statements. The easiest-to-implement hack I know of to get it to work is to use dynamic typing:

X.Foo((dynamic)Value);

Typically what I do instead is redesign the implementation so that the generic doesn't need to be used with value types. Sometimes this is pretty difficult, but often you will end up with code where the compiler is figuring everything out at compile time again, which is better.
0

##### Share on other sites

C# generics are really cumbersome to use with value types.  As far as I know you have to move the type resolution from compile time to runtime. You could make an overload of X.Foo which takes an object and then calls the appropriate method using a series of "Value is Type" if/else statements. The easiest-to-implement hack I know of to get it to work is to use dynamic typing:

X.Foo((dynamic)Value);

Typically what I do instead is redesign the implementation so that the generic doesn't need to be used with value types. Sometimes this is pretty difficult, but often you will end up with code where the compiler is figuring everything out at compile time again, which is better.

I see, that's unfortunate that c# can't do compile-time type detection for generics .  hmm, so my options are to type-cast to an object, and use reflection, or to simple create a unique instance of each type of class(which would be quite cumbersome...).

i suppose i'll use the type-cast mechanism for now, and so long as i don't see any terrible performance hit, i should be alright.

also, dynamic doesn't seem to be a known keyword for the c# compiler i'm using atm.

0

##### Share on other sites

I see, that's unfortunate that c# can't do compile-time type detection for generics .

C# can absolutely do compile time type detection for generics (specifically method parameters). That's not what you're asking it to do with this code. You're looking for C++ style duck typing and expecting the code to be generated the same way it is for c++ templates. The two are very different.

also, dynamic doesn't seem to be a known keyword for the c# compiler i'm using atm.

dynamic was introduced in .Net 4, I think. (VS2010)

Can you give a more concrete example of what you want to do? There may be better ways of solving the problem.

0

##### Share on other sites

I see, that's unfortunate that c# can't do compile-time type detection for generics .

C# can absolutely do compile time type detection for generics (specifically method parameters). That's not what you're asking it to do with this code. You're looking for C++ style duck typing and expecting the code to be generated the same way it is for c++ templates. The two are very different.

also, dynamic doesn't seem to be a known keyword for the c# compiler i'm using atm.

dynamic was introduced in .Net 4, I think. (VS2010)

Can you give a more concrete example of what you want to do? There may be better ways of solving the problem.

the code i provided is a solid example.  i'm making an interface to PSM's shaderprogram, which has ~15 overloads of various uniform inputs(int's, floats, vector's, matrixs).  i created an material class, and then i can attach nodes to that material class(each node is a uniform, or a texture). the node has a uniform location/texture index, and a value(such as an int, float, vector, or matrix), and i use the generic to pass what type the value is suppose to be.  and then when i go to apply the value to the shader, i wanted it to choose the correct overload when passing the value into the shader.

i'm open to other possibilitys, but for now i'm simply using if-else branch to check the type.

Edited by slicer4ever
0

## Create an account or sign in to comment

You need to be a member in order to leave a comment

## Create an account

Sign up for a new account in our community. It's easy!

Register a new account


text:
# Super-plethysm?

Let $U$ be a representation of $S_m$ and $V$ a representation of $S_n$. Then the representation $\operatorname{Ind}_{S_m\wr S_n}^{S_{mn}}(U^{\otimes{n}}\otimes V)$ has a nice interpretation in terms of symmetric functions. If $\operatorname{ch}(U)$ and $\operatorname{ch}(V)$ are the Frobenius characteristics of $U$ and $V$ (symmetric functions of degree $m$ and $n$), then $$\operatorname{ch}\operatorname{Ind}_{S_m\wr S_n}^{S_{mn}}(U^{\otimes{n}}\otimes V)= \operatorname{ch}(V)[\operatorname{ch}(U)],$$ where the square brackets denote plethysm. I am interested in a graded version of this result.

More specifically, let $X$ be a space with an action of $S_m$ and $Y$ a space with an action of $S_n$. I am interested in the representation $$\operatorname{Ind}_{S_m\wr S_n}^{S_{mn}}\left(H^*(X^n)\otimes H^*(Y)\right)$$ of $S_{mn}$. The naive statement would be that $$\operatorname{ch}\operatorname{Ind}_{S_m\wr S_n}^{S_{mn}}\left(H^*(X^n)\otimes H^*(Y)\right) = \operatorname{ch}(H^*(Y))[\operatorname{ch}(H^*(X))],$$ where $\operatorname{ch}(H^*(X))$ and $\operatorname{ch}(H^*(Y))$ are symmetric functions of degree $m$ and $n$ with coefficients in the polynomial ring $\mathbb{Z}[t]$.

I'm pretty sure that this statement is correct when the cohomology of $X$ is concentrated in even degree, but otherwise it is wrong. The issue is that, at the level of $S_n$-representations, the Kunneth isomorphism $H^*(X^n)\cong H^*(X)^{\otimes n}$ needs to be interpreted in terms of super-vector spaces: the action of the simple transposition $(i,i+1)$ picks up a sign when the two cohomology classes being swapped have odd degree. For example, if $X=S^1$ and $n=2$, then $H^2(S^1\times S^1)\cong H^1(S^1)\otimes H^1(S^1)$ is the sign representation of $S_2$ rather than the trivial representation.

So my question is: Is there a "super version" of plethysm for symmetric functions with coefficients in $\mathbb{Z}[t]$ for which the last displayed equation is correct?

I know that the answer is tautologically "yes"--one can just translate from symmetric functions to representations, do the induction, and translate that, and take this as the definition of super-plethysm. But I'm looking for something explicit enough to allow me to do calculations in SAGE.

• You just need to use an auxiliary variable $t$ for homological degree, and use the weight $(-t)^i$ for the $H^i$. The plethysm must act on this variable by $p_n(t)=t**n$. – F. C. Mar 16 '16 at 15:45
• Ah, that makes sense, and it indeed seems to work in a simple example. Thanks! (If you want to post that as an answer, I'll accept it.) – Nicholas Proudfoot Mar 16 '16 at 16:26
• Does this provide an answer to my MO question mathoverflow.net/questions/33681/… with t=q? – Abdelmalek Abdesselam Mar 16 '16 at 19:27

The correct way to handle these computations using plethysm is to introduce an auxiliary variable $t$ and to weight the cohomology $H^i$ with the weight $(-t)^i$. The plethysm must act on $t$ by $p_n(t)=t^n$.
• I was recently thinking about the same question. I've got a little bit different formula, but it is equivalent to the one above modulo the change of variables $q=-t$. If one introduces an auxiliary variable $q$ to weight the cohomology $H^i$ with the weight $q^i$, then the plethysm acts on $q$ by $p_n(q)=(-1)^{n-1}q^n$. – Victor May 19 '17 at 17:25

format: plain


text:
# Google wants back my microphone

My “writing” work currently goes somewhere else and have little motivation to write anything here. But, there’s something that only internet can help, whether that’s through actual help or simply transferring my annoyance.

In the past few days/weeks there has been some uproar about Facebook listening to us and later subtly suggesting products about which we talked with others. With these it’s hard to point who is objective, so I’ll paste link to web searches and I’m sure you’ll find some “evidence” – Google, Bing and DuckDuckGO. Let me also suggest Reply All podcast who recently had episode on this mysteriously called Is Facebook Spying on You?. Obviously Facebook denies all of this, but they confirm having lots of information about you whether that’s from you directly or from your friends.

Facebook and I are not in good terms for a long time. It’s more a fun social experiment rather than actual social platform. Since it isn’t on my phone there’s nothing to complain about, but there’s another omnipresent God – Google. Actually I have one of its branded phone with turned on Google Assistant, so it had to be there and had to listen to me.

Long story short, I removed microphone permissions from all Google services. Obviously some weren’t happy with this, but I can’t see how this should affect their usage. Except for Google Assistant or occasional input features, nothing should care, right? No. This is really tough break up as from time to time I’m getting vocal suggestions that are close to being commands. Google calls me to when it’s safe you’ll first need to use your phone’s screen and tap the notification then you can let the Google App access some things on your device. This is especially annoying when I’m listening to podcasts or music.

In the beginning this would go on and on, but now it’s more once a day. I don’t think that it has some “time decreasing” variable build in, so it’s definitely my action. More surprising is that even if I quickly unlock phone there won’t be anything new to give permissions to. Also, it might be only happening when the phone is locked as I haven’t had this happening otherwise.

# Free AWS is good. Not awesome, but good.

Amazon with it’s Amazon Web Service (AWS) is pretty cool. It gives you access to remote machine which you don’t have to maintain. Actually you don’t have to do anything other than use it. All machines come in different flavours, but what tastes better than free? Granted that it’s extremely limited, but surely we can squeeze something out of it. Right?

AWS instances, i.e. remote machines, differ in the amount of RAM, disc space, operating system, whether they have GPU access and so on. As you can expect free tier instance is pretty low on all measure values. To be more precise free tier instance is of t2.micro type, which is a general purpose burstable instance with a single CPU, 1 GiB memory and EBS data storage (default 4Gb storage).

What is this good for? Depending on the needs, this might be good for almost anything that doesn’t require whatever these instances are lacking. (Did I help?) Obviously. So it’s not so good for heavy computations, training machine learning models or storing data. First of all, it’s better to use for these some other services like S3, DynamoDB, Lex or general machine learning. However, in case of specific requirements, it’s always better just to rent(?) more powerful instance.

These cheap instances, in my option, are very good for few tasks. The main one is web scrapping. This is tedious task that requires small CPU bandwidth, but constant access to the internet. Moreover, we don’t really want to make many calls in small time period so there needs to be a delay between each download. That’s either because we would like to avoid being detect as a bot, or for simply politeness to the owner of the server (not clogging bandwidth).

Internet is full of examples of scrappers for different type of data. I’m adding my own to the collection with r-u-listening project. The core of the project is to allow for users to find similar music to their input. It is a bit more than recommender, but more on this project probably in the future. The scraper itself is more in two parts, i.e. crawler.py and scraper.py. The database that I’m using is FreeMusicArchive.org, which goes with slogan “It’s not just free music; it’s good music”. I do recommend it and once I have something valuable I’d like to share it with them.

Unfortunately these instances don’t come with big default memory and storage. By default they have only 4 Gb storage, which when downloading mp3 tracks will be enough for about 800 tracks (assuming about 5 Mb per track). Again, as always, it depends on the task, but for machine learning algorithms we go with The more, the merrier.

As mentioned before, free tier instances allow up to 32 Gb. To do so go to EC2 service in your AWS console. In the options tab (left side) find Elastic Block Store (EBS) and select Volumes. Then select your instance and Actions, and Modify Volume. Simple, right? In all honesty, like many things in the AWS.

I’ve been using AWS for a while. Even finished AWS general course, its essentials and 3 day onsite workshop on Architecting on AWS. All is pretty simple and consistent. I like it.

# Python Empirical Mode Decomposition on Image

One of the packages I intend long term maintain and support is Python implementation of Empirical Mode Decomposition (EMD) called PyEMD. I will skip introduction of the method as it has been explained in few other posts [1, 2, 3, …]. This blog entry is more about announcement of new feature which also means new version.

PyEMD version 0.2 is out. This means that PyEMD now supports 2D data (image) decomposition. Other visible improvements include documentation and more thorough testing both of code and data cases. Installation instructions are provided on the project’s webpage.

I am more than happy to include other improvements or suggestions. The next big step will be support for 3D and multi dimensional data. Please get in touch if you feel that there is something missing.

Image decomposition is based on the simple extremum definition: a point that is above (max) or below (min) surrounding. Behind the hood this is done using SciPy’s ndim maximum_filter. These are then connected using SmoothBivariateSpline. Stopping criteria can be chosen to be either based on the number of sifting operations or threshold values for mean and standard deviations.

Below is included exemplary decomposition, with the top image being input and the following two are the outputs. Exact formula with which the image was generated is
$sin(4\pi \cdot x) \cdot \left( cos(8\pi y + 8\pi x) + 3\right) + \left(5x + 2y - 0.4 \right) y + 2$. Python code generating this example is in provided in documentation in Examples/EMD2D.

# Easy clipboard in bash

Those who know me, know that I’m not a fan of using mouse. Whenever possible I try to avoid using it. Although, after years of experience, I am fluent in keyboarding there are still some tasks that I need mouse. Well, until quite recent.

The number one of annoying tasks is copying things from terminal to clipboard. Depending what exactly I needed to do with it I’d either select something with mouse and then copy, or stream output to file and then select it within editor. Some people found it weird, but yes, often copying through editor is much faster.

Revolution came with xclip. I’m not yet fluent in doing everything with it, but even with limited experience that I have it feels like a superpower. This program allows to copy into/from X clipboard.

Most unix distro have it installed by default or at least in package manager. In case of Ubuntu one can install it with:
$sudo apt-get update$ sudo apt-get install xclip

Examples:
Copy current directory into clipboard:
$pwd | xclip Paste whatever you have in clipboard: $ xclip -o

Copy to global clipboard so you can use in any other program:
$echo$PATH | xclip -selection clipboard

If this is too long to type, one can always use aliases, i.e. setting shorter name for long command. Either type directly into terminal (but that is only for current terminal), or update your ~/.bashrc file with:
alias "c=xclip" alias "v=xclip -o" alias "cs=xclip -selection clipboard" alias "vs=xclip -o -selection clipboard"

Then copying content of a file to a clipboard:
$cat file.txt | cs … and CTRL+V wherever one wishes. Yes, this is awesome! # Kuramoto in Python Code for Kuramoto in Python is available here or from code subpage. Explanation on how to use it is on the bottom of this post. ## Tiny introduction Kuramoto[1, 2] is probably one of the most popular and successful models for coupled oscillators. There is plenty of information about it, but in brief summary it models oscillators’ phases to be dependent on scaled phase differences for all pairs of oscillator. A while ago I have actually wrote on a specific method that I used to find optimal parameters for a Kuramoto model given some observation. (See Bayesian Dynamic Inference here.) Mathematically speaking this model is defined by a set of coupled ordinary differential equations (ODEs). Given N oscillators, dynamics for each oscillator’s phase $\phi_i$ is defined as $i\in N$ is defined by $\dot\phi_i = \omega_i + \sum_{j=0}^N k_{ij} \sin(\phi_i-\phi_j)$, where the summation is over all others oscillators. Note that one can leave coupling with itself, because in such case $\sin(\phi_i - \phi_i) = 0$, thus it doesn’t input anything into the final solution. Quick note that Kuramoto coupling model can be used to reference model with coupling function that can be represented as a sum of harmonic functions. For example, second order means that we are including both $k \sin(\Delta\phi)$ and $k^{(2)}\sin(2\Delta\phi)$. In general Kuramoto model of order $M$ would describe $\dot\phi_i = \omega_i + \sum_{m=0}^{M} \sum_{j=0}^N k_{ij}^{(m)} \sin(m (\phi_i-\phi_j))$. ## Examples Below are two examples. Both use the same core values with the exception that for the second system the coupling function defined by two harmonic terms. This naturally changes dynamics, but it shouldn’t change much average value, which is close to intrinsic frequency. Exact values for the first experiment are presented in table below. Each row is a different oscillator with initial phase Φ0 and intrinsic frequency Ω. The following columns denote coupling strength between respective pairs of oscillators. Φ0 ω k.1 k.2 k.3 1 0 28 0.2 1.1 2 π 19 0.5 -0.7 3 0 11 0.3 0.9 Phase dynamics, i.e. time derivative of obtained phases, are presented in Fig. 1. One can see that all plots are centred on intrinsic frequency with some modulations. Fig. 1. Phase dynamics, i.e. time derivative, in simple Kuramoto model. Table below shows values used in the second simulation. Extra 3 columns denote scaling values used in second harmonic. Φ0 ω k.1 k.2 k.3 k(2).1 k(2).2 k(2).3 1 0 28 0.2 1.1 -0.5 0.2 2 π 19 0.5 -0.7 -0.4 1.0 3 0 11 0.3 0.9 0.8 0.8 Again phase dynamics have been presented in a form of plot (Fig. 2). I think the difference is clear. Despite having similar mean values (approximately equal to intrinsic frequency), their modulations have change. Not only the frequency content of these modulations has changed, but also their amplitude. Fig. 2. Phase dynamics, i.e. time derivative, in Kuramoto model with included second harmonic. ## Using code Except for downloading code from either github or code subpage, one is expected to have SciPy module. Kuramoto uses it to solve differential equations. Script below shows an example of how one can use the Kuramoto module. When you run this, make sure that kuramoto.py is either in your path. Note also that most of the code is to just defining initial parameters and plotting the results. Actual execution of the module are two lines. Fig. 1 is the expected output. import numpy as np import pylab as plt from kuramoto import Kuramoto # Defining time array t0, t1, dt = 0, 40, 0.01 T = np.arange(t0, t1, dt) # Y0, W, K are initial phase, intrinsic freq and # coupling K matrix respectively Y0 = np.array([0, np.pi, 0]) W = np.array([28, 19, 11]) K1 = np.array([[ 0, 0.2, 1.1], [0.5, 0, -0.7], [0.3, 0.9, 0]]) K2 = np.array([[ 0, -0.5, 0.2], [-0.4, 0, 1.0], [ 0.8, 0.8, 0]]) K = np.dstack((K1, K2)).T # Passing parameters as a dictionary init_params = {'W':W, 'K':K, 'Y0':Y0} # Running Kuramoto model kuramoto = Kuramoto(init_params) odePhi = kuramoto.solve(T).T # Computing phase dynamics phaseDynamics = np.diff(odePhi)/dt # Plotting response nOsc = len(W) for osc in range(nOsc): plt.subplot(nOsc, 1, 1+osc) plt.plot(T[:-1], phaseDynamics[osc]) plt.ylabel("$\dot\phi_{%i}\$" %(osc+1))
plt.show()

[1] Y. Kuramoto, “Chemical Oscillations, Waves, and Turbulence,” 1984, doi: 10.1007/978-3-642-69689-3.
[2] Steven H. Strogatz, “From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators,” 2000, doi: 10.1016/S0167-2789(00)00094-4.

# Matrix Multiplication with Python 3.5

Only recently I have started to use Python 3. It’s been out for good 8+ years and all these excuses about incompatibility with some packages were just lazy. Most packages I use are already ported and if I ever find that something is incompatible… well, I’ll think then. But for now let me pat myself on the back for this great leap, because:

In Python 3.5.3 (released today) there is an operator for matrix multiplication! Check out: PEP 465 — A dedicated infix operator for matrix multiplication. The choice of operator, @, is a bit unfortunate, because of the decorators and general association with reference/internet, but seeing how few possibilities are left it’s probably the best choice.

Yes, this is big news for me. The number of times I confused myself with my own matrix operations is just too damn high! I cannot agree more with the author of the PEP 465, so let my shamelessly copy&paste (paraphrased) his reasoning. Behold!

(…) encounter many mathematical formulas that look like:

S = ( H β r ) T ( H V H T ) − 1 ( H β r )

Here the various variables are all vectors or matrices (details for the curious: [5] ).

Now we need to write code to perform this calculation. In current numpy, matrix multiplication can be performed using either the function or method call syntax. Neither provides a particularly readable translation of the formula:

import numpy as np
from numpy.linalg import inv, solve

# Using dot function:
S = np.dot((np.dot(H, beta) - r).T,
np.dot(inv(np.dot(np.dot(H, V), H.T)), np.dot(H, beta) - r))

# Using dot method:
S = (H.dot(beta) - r).T.dot(inv(H.dot(V).dot(H.T))).dot(H.dot(beta) - r)

With the @ operator, the direct translation of the above formula becomes:

S = (H @ beta - r).T @ inv(H @ V @ H.T) @ (H @ beta - r)

Notice that there is now a transparent, 1-to-1 mapping between the symbols in the original formula and the code that implements it.

Of course, an experienced programmer will probably notice that this is not the best way to compute this expression. The repeated computation of H β r should perhaps be factored out; and, expressions of the form dot(inv(A), B) should almost always be replaced by the more numerically stable solve(A, B) . When using @ , performing these two refactorings gives us:

# Version 1 (as above)
S = (H @ beta - r).T @ inv(H @ V @ H.T) @ (H @ beta - r)

# Version 2
trans_coef = H @ beta - r
S = trans_coef.T @ inv(H @ V @ H.T) @ trans_coef

# Version 3
S = trans_coef.T @ solve(H @ V @ H.T, trans_coef)

Notice that when comparing between each pair of steps, it’s very easy to see exactly what was changed. If we apply the equivalent transformations to the code using the .dot method, then the changes are much harder to read out or verify for correctness:

# Version 1 (as above)
S = (H.dot(beta) - r).T.dot(inv(H.dot(V).dot(H.T))).dot(H.dot(beta) - r)

# Version 2
trans_coef = H.dot(beta) - r
S = trans_coef.T.dot(inv(H.dot(V).dot(H.T))).dot(trans_coef)

# Version 3
S = trans_coef.T.dot(solve(H.dot(V).dot(H.T)), trans_coef)

Readability counts! The statements using @ are shorter, contain more whitespace, can be directly and easily compared both to each other and to the textbook formula, and contain only meaningful parentheses. This last point is particularly important for readability: when using function-call syntax, the required parentheses on every operation create visual clutter that makes it very difficult to parse out the overall structure of the formula by eye, even for a relatively simple formula like this one. Eyes are terrible at parsing non-regular languages. I made and caught many errors while trying to write out the ‘dot’ formulas above. I know they still contain at least one error, maybe more. (Exercise: find it. Or them.) The @ examples, by contrast, are not only correct, they’re obviously correct at a glance.

Again: yes!


text:
# LED compass

In this section, we'll implement a compass using the LEDs on the F3. Like proper compasses, our LED compass must point north somehow. It will do that by turning on one of its eight LEDs; that LED should point towards north.

Magnetic fields have both a magnitude, measured in Gauss or Teslas, and a direction. The magnetometer on the F3 measures both the magnitude and the direction of an external magnetic field but it reports back the decomposition of said field along its axes.

See below, the magnetometer has three axes associated to it.

Only the X and Y axes are shown above. The Z axis is pointing "out" of your screen.

Let's get familiar with the readings of the magnetometer by running the following starter code:

#[inline(never)]
#[no_mangle]
pub fn main() -> ! {
loop {
iprintln!("{:?}", lsm303dlhc::magnetic_field());
delay::ms(1_000);
}
}

This lsm303dlhc module provides high level API over the LSM303DLHC. Under the hood it does the same I2C routine that you implemented in the last section but it reports the X, Y and Z values in a I16x3 struct instead of a tuple.

Locate where north is at your current location. Then rotate the board such that it's aligned "towards north": its North LED should be pointing towards north.

Now run the starter code and observe the output. What X, Y and Z values do you see?

# itmdump's console
(..)
I16x3 { x: 322, y: 24, z: -11 }
I16x3 { x: 323, y: 24, z: -13 }
I16x3 { x: 324, y: 23, z: -12 }
I16x3 { x: 324, y: 20, z: -10 }

Now rotate the board 180 degrees. What X, Y and Z values do you see this time?

Finally, rotate the board 90 degrees clockwise. The East LED should be pointing north this time. What X, Y and Z values do you see this time?


text:
# Talk:Autoregressive model

WikiProject Statistics (Rated C-class, High-importance)

This article is within the scope of the WikiProject Statistics, a collaborative effort to improve the coverage of statistics on Wikipedia. If you would like to participate, please visit the project page or join the discussion.

C  This article has been rated as C-Class on the quality scale.
High  This article has been rated as High-importance on the importance scale.

## Summary and Better Examples?

This article is noticeably more difficult to comprehend than many other Wiki entries on statistics. A higher level summary statement on AutoRegressive techniques could be useful, especially example usages that point the reader in the right direction (ARMA, ARCH, ....) .

Strawman example of opening statement: Autoregressive models are used for prediction and data smoothing, especially in time-series data containing signal noise(SNR), moving averages(ARMA), and characteristically different time periods(ARCH). AutoRegressive (AR) models are linear combinations of the input parameter values and white noise values.

As an example, consider an input dataset with two parameters, time and power. We want to predict how much power we need to provide in the next second. Some noise exists in the power measurement, but we dont know exactly how much. We also know that real power usage varies considerably, for example during peak usage periods vs off-peak energy usages. Using an autoregressive model, we can learn from previous instances and predict the next power value. (Picture of example) — Preceding unsigned comment added by 208.127.244.182 (talk) 00:22, 30 November 2012 (UTC)

## Constant term

Some of the formulas include the constant term:

$X_t = c + \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t.\,$
while other formulas don't:
$X_t = \sum_{i=1}^p \varphi_i X_{t-i}+ \varepsilon_t.\,$
I think the article should be consistent in the notation. Albmont (talk) 18:26, 18 November 2008 (UTC)
You're absolutely right. Personally I have never encountered the constant term outside of Wikipedia, but that's just me. --Zvika (talk) 06:43, 19 November 2008 (UTC)
gretl includes the constant term (optionally). So does R (programming language), as in fit.ar.par. But both softwares include anything in the "deterministic" component. Albmont (talk) 12:21, 19 November 2008 (UTC)
Well, then, definitely go ahead and put it in. --Zvika (talk) 13:39, 19 November 2008 (UTC)

I think a problem with the constant term is that the equations for wide-sense stationarity involving the poles of the time-shift polynomial are no longer applicable. Consider the AR(1) model with the constant $c$ term. Assuming some nonzero initialization for the sequence, it's trivial to expand the sequence at a given time and use a geometric series to show that

$\mathbb{E}[X_t] = c \frac{1-\phi^n}{1-\phi} + \phi^{n-1} X_0,$

which varies with $n$ but converges to $c/(1-\phi)$ in the limit. In other words, for the system to be wide-sense stationary, $c$ must equal zero. This result is compatible with Exercise 1.6 in "Adaptive Filter Theory" (4th Ed.) by Simon Haykin, which states that the input to an AR(1) must have zero mean. I believe this is an error in both this article and the general ARMA article. — Preceding unsigned comment added by 152.3.43.164 (talk) 18:02, 11 March 2013 (UTC)

## Autocovariance or autocorrelation?

According to http://en.wikipedia.org/wiki/Spectral_density the spectral density is the FT of the autocorrelation (and according to my notes!) but here it is stated that it is the FT of the autocovariance. In the case there μ = 0 it doesn't affect the result, but is it right? If so can someone clarify the apparently contradicting information? —Preceding unsigned comment added by 163.1.167.139 (talk) 22:43, 21 March 2009 (UTC)

## Gretl

Apropos gretl, I just noticed that when gretl computes the parameters in the AR(1) model with a constant term, it returns const and phi_1 based on equation $X_t = c + \varphi (X_{t-1} - c) + \epsilon_t\,$ instead of $X_t = c + \varphi X_{t-1} + \epsilon_t\,$. Albmont (talk) 17:05, 19 May 2009 (UTC)

## Variance of X_t

The variance of Xt should be only valid for the assymptotic case, when t goes to infinity. That value is valid for a finite t in Xt only when we have a process that begins at t = minus infinite; in most real-world applications (for example, Monte Carlo simulations of AR(1) series), we begin with a fixed X0 and, depending on phi and sigma, we may never get even close to the assymptotic values. Albmont (talk) 13:47, 9 October 2009 (UTC)

At least in the textbooks that I use, an AR process is defined as one which begins at negative infinity (e.g., Porat's "Digital Processing of Random Signals"). This is required to ensure that the process is wide-sense stationary. Almost all of the text of the article would change if you were to change this definition. For example, it would no longer be possible to talk about the autocovariance of the process or its spectral density. --Zvika (talk) 14:20, 9 October 2009 (UTC)
Maybe it could be possible to get a compromise. Let's write non-assymptotic equations for the conditional AR(1), namely Xt|X0 - or even Xt|Xs, t > s. I think these formulas are more useful (for the sake of Monte Carlo analysis) than the assymptotic equations for a hypothetical series that begins at t = -infinite. Albmont (talk) 14:26, 9 October 2009 (UTC)
I don't object to that in principle, if it is stated in addition to the existing formulas, say in a separate section on conditional properties. --Zvika (talk) 15:25, 9 October 2009 (UTC)
OTOW, the current version of the article implies some inconsistencies. Because it says that AR(1) = random walk for φ = 1 (this is true only when X0 is zero), also it allows an AR(p) even when the coefficients have a unit root (or worse). Maybe it should be better to keep the analysis of the stationary process with |φ| < 1 and t beginning at -infinite at a separate section too. Albmont (talk) 16:07, 9 October 2009 (UTC)

(outdent) OK, I was not aware of the fact that a random walk necessarily equals 0 at time 0, but apparently this is what it says in the random walk article, so I reworded that part. I still maintain that the standard definition of an AR process begins at -infinity. Do you have a source that says something else? --Zvika (talk) 08:55, 10 October 2009 (UTC)

## State space form

AR(p) model $X_{t}=\phi_{1}X_{t-1}+\phi_{2}X_{t-2}+\cdots+\phi_{p}X_{t-p}+\varepsilon_{t},\; t\ge p$ where $\varepsilon_{t}\sim N\left(0,\sigma^{2}\right)$

$X_{t}=X_{t},\; t

The usual estimation method doesn't full use the data points t=0,...,p-1. Introduce the state space form, and to some extent, we can use these data points in a better way.

define $e_{1}^{'}=\left(\begin{array}{cccc} 1 & 0 & \cdots & 0\end{array}\right), E_{t}^{'}=\left(\begin{array}{cccc} X_{t} & X_{t-1} & \cdots & X_{t-p+1}\end{array}\right)$

$G=\left(\begin{array}{ccccc} \phi_{1} & \phi_{2} & \cdots & \phi_{p-1} & \phi_{p}\\ 1 & 0 & \cdots & 0 & 0\\ 0 & 1 & \cdots & 0 & 0\\ 0 & 0 & \ddots & 0 & 0\\ \vdots & \vdots & \cdots & \vdots & \vdots\\ 0 & 0 & \cdots & 1 & 0\end{array}\right)=\left(\begin{array}{cc} \begin{array}{cccc} \phi_{1} & \phi_{2} & \cdots & \phi_{p-1}\end{array} & \phi_{p}\\ I_{p-1} & 0\end{array}\right)$,

then the state space form is

$X_{t}=e_{1}^{'}E_{t},\; t\ge p$

$\boldsymbol{\boldsymbol{X}}_{p}=G\boldsymbol{\boldsymbol{X}}_{p-1}^{'}+e_{1}\varepsilon_{p} where \boldsymbol{\boldsymbol{X}}_{p}^{'}=\left(\begin{array}{cccc} X_{p-1} & X_{p-2} & \cdots & X_{0}\end{array}\right)$

$\text{E}\left(\boldsymbol{\boldsymbol{X}}_{p}\right)=0$ $\text{Var}\left(\boldsymbol{\boldsymbol{X}}_{p}\right)=G\text{Var}\left(\boldsymbol{\boldsymbol{X}}_{p-1}\right)G^{'}+\text{Var}\left(\varepsilon_{p}\right)e_{1}e_{1}^{'}$, if stationarity is imposed, then $\text{Var}\left(\boldsymbol{\boldsymbol{X}}_{p}\right)=\text{Var}\left(\boldsymbol{\boldsymbol{X}}_{p-1}\right)=\Omega_{p}$, i.e.$\Omega_{p}=G\Omega_{p}G^{'}+\sigma^{2}e_{1}e_{1}^{'}.\text{vec}\left(\Omega_{p}\right)=\left(G\otimes G\right)\text{vec}\left(\Omega_{p}\right)+\sigma^{2}\text{vec}\left(e_{1}e_{1}^{'}\right)=\sigma^{2}\left(I-G\otimes G\right)^{-1}\text{vec}\left(e_{1}e_{1}^{'}\right).$

$\text{vec}\left(ABC\right)=\left(C^{'}\otimes A\right)\text{vec}\left(B\right).$ Jackzhp (talk) 19:25, 26 March 2011 (UTC)

## OLS procedure

I feel that it is necessary to mention the reason why people don't apply the Ordinary least squares to estimate the coefficients. Jackzhp (talk) 19:25, 26 March 2011 (UTC)

## AR(2) Spectrum

The page used to state that

For AR(2), the spectrum has a minimum ($\varphi_2 > 0$) or maximum ($\varphi_2<0$) if[citation needed]
$|\varphi_1(1-\varphi_2)| < 4|\varphi_2|.$

However, I am almost certain this is wrong. The critical points of the AR(2) spectrum occur when

$\varphi_1(1-\varphi_2)\sin(\omega) + 4\varphi_2\sin(\omega)\cos(\omega) = 0$

Thus they occur at $\omega = 2\pi k$ or at $\omega = \cos^{-1}\frac{\varphi_1(1-\varphi_2)}{4\varphi_2}$. I believe the person who posted the above made the mistake of dividing by sin (which is sometimes zero) and thus eliminating some of the potential peaks.

AR(0), AR(1), and AR(2) processes with white noise

This graph caption in the section "Graphs of AR(p) processes" is inadequate. There are five subgraphs but the caption only explains three (presumably the top three??). The numbers by the right side of each subgraph are undefined. And from the original documentation, I can't even confirm that the caption is correct in referring to the top three subgraphs. Can someone figure this out and redo the caption? Thanks. Duoduoduo (talk) 15:31, 10 January 2013 (UTC)

AR(0), AR(1), and AR(2) processes with white noise

The graph makes sense in the context of the section where it is included. There is only one possible plot for AR(0) since there are no parameters. There are two plots for AR(1), one for a value of φ close to zero and another for φ just less than one. The last two plots are for AR(2). One plot is for where φ1 and φ2 have the same sign. The other plot is when the two parameters have different signs.

I agree the graph could use some more work to make it clearer. I will look into it.— Preceding unsigned comment added by Everettr2 (talkcontribs) 01:23, 13 January 2013 (UTC)

Thanks. I've clarified the caption based on your explanation. Duoduoduo (talk) 01:42, 13 January 2013 (UTC)

## Edit to lede

@TheSeven: Please keep in mind WP:BRD -- you boldly made an edit, deleting a passage; I reverted your edit, then you are supposed to discuss on the talk page rather than edit war.

Your edit in the two-sentence long lede changes the second and last sentence from

The autoregressive model is one of a group of linear prediction formulas that attempt to predict an output of a system based on the previous outputs.

to

The autoregressive model is one of a group of linear prediction formulas.

thereby removing the passage

that attempt to predict an output of a system based on the previous outputs.

the wording suggested that they are only for prediction, which is untrue; also, prediction is already mentioned--keep intro clear and simple

(1) Not sure what you have in mind about other uses of AR. It seems to me that any others must be very minor or must actually be aspects of prediction. Can you be specific?

(2) I agree that predict should not be in there twice.

(3) The first sentence of the lede says that the AR model is a process, while the second sentence (both with and without your edit) says that it's a formula. That's awkward and needs to be fixed.

(4) Your deletion of the passage removes from the lede the most important thing there is to say about AR: based on the previous outputs. We can't possibly have a lede that doesn't even mention that.

I'm going to revert your edit as a violation of BRD and then rewrite the lede to take these things into account. Feel free to discuss here or to tweak my new version, or even to revert my new version to the original version, but please don't restore your version unless and until there arises a consensus on the talk page to do so. Duoduoduo (talk) 17:39, 30 January 2013 (UTC)

I really like the most-recent version. TheSeven (talk) 17:44, 31 January 2013 (UTC)


text:
# Thread: find the value for which the function is continuous

1. ## find the value for which the function is continuous

f(x)= ln x/(x-1), if 0<x not =1
c , if x=1
and what kind of discontinuity is present if c does not have this value?

can someone help me with this problem?

2. This is nearly impossible to read...

Is it

$f(x) = \begin{cases}\frac{\ln{x}}{x - 1}\textrm{ if }x \in (0, 1)\cup (1, \infty)\\c\textrm{ if }x=1\end{cases}$

or

$f(x) = \begin{cases}\ln{\left(\frac{x}{x - 1}\right)}\textrm{ if }x \in (0, 1)\cup (1, \infty)\\c\textrm{ if }x=1\end{cases}$?

3. $f(x) = \begin{cases}\frac{\ln{x}}{x - 1}\textrm{ if } 0

thats how it is i dont know how to do a not equal sign

4. Originally Posted by dat1611
$f(x) = \begin{cases}\frac{\ln{x}}{x - 1}\textrm{ if } x>0 \, ; \, x \ne 1 \\c\textrm{ if }x=1\end{cases}$
fify

5. Originally Posted by dat1611
$f(x) = \begin{cases}\frac{\ln{x}}{x - 1}\textrm{ if } 0

thats how it is i dont know how to do a not equal sign
It is done by \neq

[LaTeX ERROR: Convert failed]

For the question at hand, you need to determine a value c such that the function is continuous. One way could be to just plug in 1 in the above definition with the natural log, but you cannot do that since you'd get something of the form 0/0. Hence you need to find the limit of [LaTeX ERROR: Convert failed] as x approaches 1. Or more succinctly, find

[LaTeX ERROR: Convert failed]

Once you find that, you can just set c equal to it. From the definitions of continuous functions, f(x) is continuous at point a if

[LaTeX ERROR: Convert failed]

6. ok can someone help with the actual problem

7. Do you know L'Hopital's rule? If you do, just apply to the numerator and denominator and it ' should be pretty easy.

8. Hello, dat1611!

$f(x)\;=\; \left\{ \begin{array}{cc}\dfrac{\ln x}{x-1}& x > 0,\;x \ne 1 \\ \\[-3mm]
c & x=1 \end{array} \right$

Find the value of $c$ for which the function is continuous.

To be continuous at $x = 1$, we want: . $\displaystyle{\lim_{x\to1}\frac{\ln x}{x-1} \;=\;c }$

The left side goes to $\frac{0}{0}$ so we can apply L'Hopital.

. $\displaystyle{\lim_{x\to1}\frac{\frac{1}{x}}{1} \;=\;c \quad\Rightarrow\quad \lim_{x\to1}\frac{1}{x} \:=\:c$

Therefore: . $c \:=\:1$

9. so
$\displaystyle \lim\limit_{x \to 1}\frac{\ln{x}}{x - 1}
$
= 1


text:
Intel Parallel Studio for improving OpenCV

I would like to contribute to OpenCV making a parallel version (using OpenMP) of SURF.

Yes, yes, I know, SURF in OpenCV is already "parallel" and it uses Intel TBB, but come on it scales horribly, I'm sure some improvement can be done for sure ;)

I want to use Intel Parallel Studio tools (e.g. VTune Amplifier, Intel Advisor and Intel Inspector) for making this process not a nightmare and make the paralelization easier. However, for some of them (Intel Advisor for sure) we need to include some libraries in the src code (e.g. #include <advisor-annotate.h>) and of course editing the makefile. For example, this is the compile command for one of the Intel Advisor samples using in the tutorial:

g++ -m64 nqueens_serial.cpp -o 1_nqueens_serial -O2 -g -I /opt/intel/advisor_2017.1.1.486553/include/ -ldl

Now there is my question (and I'm afraid that the answer will be painful): OpenCV uses CMAKE and I never developed something using it (only used it for building) and I have no idea how could I integrate some Intel includes and flags necessary for this kind of profiling/improvements in the CMAKE files or if this is even possible.

Someone can help me or give me some guidelines?

edit retag close merge delete

... a wasps nest just like you say. Without a decent knowledge of the inner CMAKE pipeline this will be hell ...

( 2017-02-03 08:52:41 -0500 )edit

So happy to hear the opencv is so easy to contribute (lol). I think that "start googling about learning CMAKE" is the best suggestion that you can give me @StevenPuttemans , right? Or do you have anything better? :)

( 2017-02-03 09:02:33 -0500 )edit

You are digging right into one of the most complex tasks in the OpenCV code i guess... There is a reason why i try to leave the cmake as is. But I guess learning CMAKE basics is indeed a good start.

( 2017-02-03 13:29:56 -0500 )edit


text:
# Convergence in mean square from expected value/variance

I'm looking for a proof of the following statement: Given a sequence of independent random variables $X_n$ satisfying $$\lim_{n\to \infty} E[X_n] = T,$$ where T is a constant, then
$$\lim_{n\to \infty} V[X_n] = 0$$ implies convergence of $X_n$ to $T$ in the mean-square. This statement is supplied without proof or reference in Shreve's Stochastic Calculus book.

-
@t-laarhoven: "Mean-square" convergence means $L^2$ convergence, i.e. we want to show $E[(X_n - T)^2] \to 0$. – Nate Eldredge Aug 24 '11 at 14:59

I assume $V[X_n]$ is the variance.
Let $\mu_n = E[X_n]$ for convenience, and write \begin{align*} E[(X_n - T)^2] &= E[(X_n - \mu_n + \mu_n - T)^2] \\ &= E[(X_n - \mu_n)^2] + (\mu_n - T)^2.\end{align*}
(The cross term vanished since $E[X_n - \mu_n]=0$.) Now both terms go to 0 by assumption.
Hmm, this is very simple. I think my mind was rebelling at the implication that this implies almost-sure convergence for a subsequence of ${X_n}$. The latter fact seems too good to be true... – James Davidoff Aug 24 '11 at 15:16
By the way, this is why $E(X)$ is the value of $x$ which minimizes $E((X-x)^2)$ (and why, as a consequence, the minimal value is the variance of $X$). – Did Aug 28 '11 at 20:45


text:
The MathNet Korea
Information Center for Mathematical Science

### 논문검색

Information Center for Mathematical Science

#### 논문검색

Tohoku Mathematical Journal. Second Series
( Vol. 58 NO.3 / (2006))
Equivariant completions of toric contraction morphisms

Pages. 303-321
Abstract We treat equivariant completions of toric contraction morphisms as an application of the toric Mori theory. For this purpose, we generalize the toric Mori theory for non-$mathcal Q$-factorial toric varieties. So, our theory seems to be quite different from Reid's original combinatorial toric Mori theory. We also explain various examples of non-$mathcal Q$-factorial contractions, which imply that the $mathcal Q$-factoriality plays an important role in the Minimal Model Program. Thus, this paper completes the foundation of the toric Mori theory and shows us a new aspect of the Minimal Model Program. 1. Introduction 2. Equivariant completions of toric contraction morphisms 3. Applications of equivariant completions 4. Examples of non-$\mathcal Q$-factorial contractions Toric varieties, Mori theory, minimal model program, equivariant completion 14M25, 14E30


text:
This article introduces predict(), fitted(), residuals() for in-sample and out-of-sample data. I will also show how to get creative with mcp, including how to make predictions around future change points.

# Preparation: an example model

We need an mcpfit to get started. We take the “demo” dataset:

library(mcp)
options(mc.cores = 3)  # Speed up sampling
set.seed(42)  # Make the script deterministic

ex = mcp_example("demo")
head(ex$data) ## # A tibble: 6 x 2 ## time response ## <dbl> <dbl> ## 1 68.4 32.8 ## 2 87.3 -1.16 ## 3 69.0 27.6 ## 4 11.6 10.1 ## 5 19.5 14.1 ## 6 46.1 18.3 … and model it as three segments, i.e., two change points: # Define the model model = list( response ~ 1, # plateau (int_1) ~ 0 + time, # joined slope (time_2) at cp_1 ~ 1 + time # disjoined slope (int_3, time_3) at cp_2 ) # Fit it fit = mcp(model, data = ex$data)

This is what the data and the inferred fit looks like with 95% credible interval and a 80% prediction interval:

plot(fit, q_fit = TRUE, q_predict = c(0.1, 0.9))

To review what we see here:

• The black dots is the data from ex$data. • The gray lines are 25 samples from the posterior (control using plot(fit, lines = 100)). • The dashed red lines are the 2.5% and 97.% quantiles of the fitted (expected) values. • The green lines are the 10% and 90% quantiles of the predicted values. • The blue curves on the x-axis are the posterior distributions of the change point locations (better viewed using plot_pars(fit, pars = c("cp_1", "cp_2"))). Behind the scenes, plot() merely calls predict() and fitted() to show these inferences. # Extracting fitted() values To get the fitted values for each data point, simply do fitted(fit): head(fitted(fit)) ## time fitted error Q2.5 Q97.5 ## 1 68.35820 30.352795 1.0524240 28.281702 32.375971 ## 2 87.29038 -3.316264 0.7639949 -4.833468 -1.826587 ## 3 69.01173 30.704731 1.0839745 28.566277 32.786091 ## 4 11.59361 10.302963 0.7323588 8.826560 11.720457 ## 5 19.50091 10.311477 0.7238347 8.878012 11.722323 ## 6 46.12009 18.377166 1.0029739 16.038041 20.057603 In general, this output will include: • A column for each predictor column in the data. Here, time is the only predictor and you see the values in the same order as in ex$data (which is copied to fit$data). But models with varying change points will additionally have a column for the varying data, binomial() models include the number of trials, etc. • fitted: The fitted value which is the mean of the posterior for this predictor value (or set of predictor values). I.e., it is the Expected Value (EV). • error:: The standard error corresponding to fitted, i.e., diff(fitted + c(-1, 1) * error) is the 68% credible interval. • Q[some number]: The quantiles of the fitted distribution. You can set the quantiles using fitted(fit, probs = c(0.1, 0.5, 0.9)). If you compare these values to the plot, you will see that they correspond. plot() merely calls fitted() behind the scenes. ## Expected values for out-of-sample data To predict out-of-sample data, you can simply use the newdata argument. newdata = data.frame(time = c(ex$data$time[1], 25, -20, 200)) fitted(fit, newdata = newdata) ## time fitted error Q2.5 Q97.5 ## 1 68.3582 30.35279 1.0524240 28.281702 32.375971 ## 2 25.0000 10.34998 0.7175872 8.967093 11.770927 ## 3 -20.0000 10.30119 0.7382389 8.825298 11.720457 ## 4 200.0000 -28.14998 10.6442779 -49.429380 -7.991097 Note that: • We get one row per value of time. • The first value for time is in the dataset. The values correspond to the same row in fitted(fit) because that’s merely a shortcut to do fitted(fit, newdata = fit$data).
• The second value (time = 20) is within the observed region, but not in the dataset.
• The third value (time = - 20) is outside the observed region, but mcp merely extends the first segment backwards in time. Because it’s a plateau, we see approximately the same values as for time = 20.
• The fourth value (time = 200) is way outside the observed region. Because it is the extrapolation of the slope in the third segment of which we’ve only observed the first tiny bit, the posterior distribution is very wide because even a small uncertainty in the slope results in very large differences further out.

## Arguments

If you look at the documentation for predict(), fitted(), and residuals(), you’ll see that they are quite versatile, taking many different arguments. To mention a few, you can set which_y = "sigma" to get fitted values for sigma more on modeling sigma, prior = TRUE to predict using only the prior, varying = TRUE and arma = FALSE can be toggled to include/exclude AR(N) and varying effects.

# Predictions

predict() is the posterior predictive and it takes exactly the same arguments as fitted(). This means that you can make predictions for in-sample and out-of-sample data as well. As with fitted(), plot() uses predict() under the hood to plot prediction intervals. You can see that the values correspond to dashed green lines in the plot (the 80% prediction interval):

set.seed(42)
head(predict(fit, probs = c(0.1, 0.9)))
##       time   predict    error       Q10       Q90
## 1 68.35820 30.275914 4.184053 24.972560 35.608088
## 2 87.29038 -3.339362 4.101321 -8.619435  1.880431
## 3 69.01173 30.764242 4.148573 25.449155 36.020085
## 4 11.59361 10.269686 4.077599  5.100657 15.494820
## 5 19.50091 10.195011 4.107429  4.956546 15.454846
## 6 46.12009 18.388889 4.124531 13.147433 23.702492

Note that predict() uses random sampling under the hood, so these values will differ slightly from call to call. You can make it replicable using set.seed() as above. In general, the more posterior samples, the less the call-to-call variance will be. Conversely, fewer samples means more call-to-call variation, e.g., if you do predict(fit, nsamples = 10)).

# Residuals

## Step 2: add the unobserved segment(s) and fit

Now we extend the model with the future segment of which we have prior knowledge:

model_forecast = c(fit$model, list( ~ 1 # intercept (int_4) after cp_3 )) And finally, we extend the list of priors with the two new parameters (time_4 and cp_3). It may be helpful to review the article on priors in mcp. prior_forecast = c(fit$prior, list(
))
fit_forecast = mcp(model_forecast, data = ex$data, prior = prior_forecast) ## Step 3: predict! We can go right ahead and compute our 50% and 80% prediction intervals at time = 125: predict(fit_forecast, newdata = data.frame(time = 125), probs = c(0.1, 0.25, 0.75, 0.9)) ## time predict error Q10 Q25 Q75 Q90 ## 1 125 3.946098 10.99492 -13.8512 -5.578755 11.79398 14.61287 To really understand what’s going on here, it may be helpful to visualize the model. For now, we will have to hack this a bit too, manually doing our plot: # Get posterior and posterior predictive "predictions" newdata = data.frame(time = 1:170) fitted_forecast = fitted(fit_forecast, newdata = newdata, summary = FALSE, nsamples = 50) predict_forecast = predict(fit_forecast, newdata = newdata, summary = FALSE) # Plot it library(ggplot2) ggplot(predict_forecast, aes(x = time, y = predict)) + # Prediction intervals and line at x = 125 stat_summary(fun.data = median_hilow, fun.args = list(conf.int = 0.8), geom = "ribbon", alpha = 0.2) + stat_summary(fun.data = median_hilow, fun.args = list(conf.int = 0.5), geom = "ribbon", alpha = 0.3) + geom_vline(xintercept = 125, lty = 2, lwd = 1) + # Lines for fitted draws geom_line(aes(y = fitted, group = .draw), data = fitted_forecast, alpha = 0.2) + # Observed data geom_point(aes(x = time, y = response), data = ex$data) +
labs(title = "Predicting with future change points")
You can read the predicted values from above at x = 125 off this graph. We literally just predicted for all values between 1 and 170, and visualized it using a ribbon. This means that you can also predict further into the future, if you’d like.
Without doing this formal model of the future change point, one may have thought that the change point should occur around time = 110 since that’s the expected value of cp_2 + (cp_2 - cp_1). However, we truncated the prior for the future change point (cp_3) so that it occurs after the last data point (MAXX), i.e., at time > 100. This is knowledge that the third change point had not yet been observed at time = 100, and this pushes the distribution further into the future (actually around 118; see summary(fit_forecast)).


text:
## Precalculus (6th Edition) Blitzer

Identity for $\tan \left( \theta -\phi \right)$ is expressed as $\frac{\tan \theta -\tan \phi }{1+\tan \theta \tan \phi }$.
From the difference formula of tangents, the tangent of the difference between two angles, say A and B, is expressed as, $\tan \left( A-B \right)=\frac{\tan A-\tan B}{1+\tan A\tan B}$ Thus, for angles $\theta$ and $\phi$ , $\tan \left( \theta -\phi \right)=\frac{\tan \theta -\tan \phi }{1+\tan \theta \tan \phi }$